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1-1.  Oxidation/Reduction  Behavior  of  Pure  Indium 

Overview 

This  work  focused  on  the  study  of  oxidation  and  reduction  of  pure  indium  solder 
to  subsequently  develop  a  reflow  processing  window  to  be  utilized  in  industry.  Work 
began  with  a  thermodynamic  survey  in  which  the  Gibbs  free  energy  of  oxidation  was 
calculated.  From  there,  the  effect  of  moisture  and  oxygen  on  the  hydrogen  reducing 
environment  was  investigated  to  develop  an  ideal  processing  environment. 

Before  beginning  the  oxidation  experiments  it  was  necessary  to  develop  a 
polishing  technique  which  would  yield  repeatable  thickness  measurements  for  ultra  thin 
oxide  layers.  An  electrochemical  polishing  approach  was  adopted  because  of  its  ability 
to  produce  an  ultra  smooth  (<  10-nm  roughness)  finish.  Spectroscopic  ellipsometry  was 
employed  to  measure  the  thickness  of  these  ultra  thin  oxide  layers.  Ellipsometry  is 
preferred  because  of  its  high  resolution  (5  angstrom)  and  non-destructive  nature.  For 
thicker  oxide  layers,  measurements  were  estimated  using  dynamic  nanoindentation. 

To  experimentally  verify  the  theoretical  thermodynamic  calculations,  oxide 
thickness  verses  temperature  plots  were  constructed.  From  these  plots  several 
conclusions  were  drawn.  First,  there  is  little  or  no  oxidation  when  indium  is  heat  treated 
below  the  melting  temperature.  Second,  based  on  the  Arrhenius  plot  there  is  a  sharp 
increase  in  the  activation  energy  after  crossing  the  melting  temperature.  Third,  oxidation 
in  inert  environments  seems  to  follow  thermodynamic  calculations.  A  cross  over  was 
seen  from  oxidation  to  reduction. 

Reduction  of  thick  indium  oxide  layers  in  a  hydrogen  environment  seemed  to  be 
more  difficult.  Reduction  environmental  conditions  close  to  thermodynamic  equilibrium 
point  showed  very  slow  kinetics.  Because  indium  oxide  is  so  stable  at  high  temperatures 
it  was  necessary  to  use  extreme  temperatures  (350°C)  to  see  any  appreciable  reduction. 

From  a  kinetics  standpoint  it  was  seen  that  oxidation  growth  in  air  follows  a 
logarithmic  relationship.  Indium  solder  heat  treated  at  145°C,  180°C,  220°C  all  showed 
this  behavior  which  is  based  on  the  theory  of  electron  flow  from  the  metal  to  the  oxide 
rather  than  diffusion  of  ions  through  the  oxide  lattice.  Oxidation  kinetics  in  an  inert 
environment  below  the  melting  temperature  showed  little  or  no  oxide  growth,  all  falling 
within  the  range  of  the  native  oxide  layer. 

Complementing  this  oxidation  and  reduction  study  is  a  materials  characterization 
of  both  indium  and  indium  oxide.  The  mechanical  properties,  elastic  modulus  and 
hardness  were  measured  using  nanoindentation.  Microstructure  characterization  was 
carried  out  using  both  optical  microscopy  and  atomic  force  microscopy  (AFM).  Lastly, 


X-ray  diffraction  (XRD)  was  used  to  verify  the  crystal  structure  of  indium  oxide  before 
and  after  the  melting  temperature. 

Thermodynamics 

Oxidation  begins  when  AG°  of  the  reaction  takes  on  a  negative  value  due  to  the 
metal  being  unstable  relative  to  the  oxide  formed  at  elevated  temperatures.  In  Fig.  1-1 
below  it  can  be  clearly  seen  that  indium  oxide  is  stable  from  0°K  to  indium's  boiling 
point  (2080°C)  meaning  that  indium  oxide  is  always  oxidizing  at  all  relevant 
temperatures.  The  only  way  to  make  the  oxide  thermodynamically  unstable  is  to 
introduce  a  reducing  gas  such  as  CO2,  or  H2. 

To  achieve  thermodynamic  equilibrium  hydrogen  gas  is  added  to  the  environment 
until  the  ratio  of  moisture  to  hydrogen  at  a  specific  temperature  is  small  enough  to  make 
AG°  =  0.  Once  AG°  is  zero  more  hydrogen  is  added  to  enter  the  reduction  zone  where  the 
oxide  is  unstable  relative  to  the  metal.  So  it  is  immediately  obvious  that  F^O  partial 
pressure  plays  a  large  role  in  determining  whether  an  environment  is  reducing  or 
oxidizing.  Moisture  should  be  reduced  as  much  as  possible  and  hydrogen  increased  to 
maximize  the  etching  rate  in  reducing  environments.  It  should  be  noted  here  that 
thermodynamics  will  only  tell  us  whether  the  environment  is  reducing  or  oxidizing.  No 
growth  or  etching  rate  information  can  be  obtained  from  such  a  study. 
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Fig.  1-1  Ellingham  diagram  for  indium  oxide,  Gibbs  free  energy  verses  temperature. 

The  graph  below  (Fig.  1-2)  of  percent  hydrogen  verses  temperature  with  moisture 
isobars,  illustrates  the  effect  of  moisture  on  the  thermodynamic  location  of  the  reduction 
zone.  We  can  see  that  as  the  moisture  content  increases  the  reduction  zone  shifts  further 
and  further  to  the  right.  Therefore  at  higher  moisture  concentrations  higher  temperatures 
are  needed  to  achieve  reduction.  For  example  at  200°C  and  4%  hydrogen,  only 
environments  with  <  lppm  moisture  are  reducing.  With  lOppm  moisture  and  4% 


hydrogen  the  temperature  would  have  to  be  increased  to  300°C  to  enter  the  reduction 
zone.  Since  processing  temperatures  greater  than  300°C  are  impractical  because  of 
component  degradation,  low  moisture  concentrations  are  imperative  in  the  fluxless  reflow 
of  pure  indium  solder. 


Fig.  1-2  Effect  of  moisture  on  the  hydrogen  reducing  environment. 


Fig.  1-3  Equilibrium  moisture  hydrogen  ratio  verses  temperature  with  oxygen  isobars. 

In  contrast  the  effect  of  oxygen  partial  pressure  doesn’t  have  much  of  an  impact 
on  the  thermodynamic  location  of  the  oxidation  and  reduction  zones.  In  Fig.  1-3  it  can  be 
seen  that  since  the  equilibrium  oxygen  partial  pressure  is  so  low  (which  correlates  to  a 
very  small  moisture  to  hydrogen  ratio)  changes  in  atmospheric  oxygen  partial  pressures 


in  the  range  of  1  O'25  atm  to  air  have  little  effect  on  the  overall  moisture  to  hydrogen  ratio. 
This  is  due  to  the  logarithmic  relationship  between  the  moisture  and  hydrogen  ratio  and 
temperature. 

Kinetics 

To  this  point  we’ve  only  discussed  whether  a  reaction  is  oxidizing  or  reducing, 
now  we  must  discuss  the  rate  at  which  these  reactions  occur.  There  are  several  kinetics 
models  for  oxidation  which  functionally  describe  how  the  oxides  grow.  These  models 
can  be  narrowed  down  to  three  main  categories,  parabolic,  inverse  logarithmic,  and 
logarithmic. 

Parabolic  growth  rate  kinetics  can  be  most  easily  described  as  when  the  diffusion 
of  ions  through  oxide  is  the  rate  controlling  mechanism.  Wagner’s  oxidation  equation 
below  shows  the  oxides  growth  dependence  on  time. 

d2  =  Ale-Q,kTt  +  A2  (i) 

where  d  is  oxide  thickness,  Q  is  activation  energy,  k  is  a  rate  constant,  t  is  time  and  A| 
and  A2  are  constants.  Cabrera  and  Mott’s  theory  of  oxidation  says  that  the  growth  of  an 
oxide  layer  follows  an  inverse  logarithmic  relationship.  Here  the  rate  controlling 
mechanism  is  the  flow  of  metal  ions.  Growth  of  the  native  oxide  layer  on  many  metals 
such  as  copper,  iron  and  indium  can  be  described  by  this  mechanism. 

Finally  there  is  Uhlig’s  theory  electron  transport  which  is  governed  by  a 
logarithmic  growth  rate.  Figure  1-4  shows  the  initial  stages  of  oxidation  where  first 
oxygen  is  physically  adsorbed  on  the  metal  surface,  attached  only  by  weak  van  der  Waals 
bonds.  Subsequently  the  oxygen  molecules  are  chemisorbed  with  the  oxygen  molecules 
acting  as  the  electron  acceptor  and  the  metal  acting  as  the  electron  donor.  Lastly  that 
metal-oxide  complex  is  sublimated  to  form  the  oxide  lattice. 
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Fig.  1-4  Illustration  of  oxidation  process. 


Here  the  rate  controlling  mechanism  is  the  follow  of  electrons  from  the  metal  to 
the  oxide.  When  enough  electrons  have  passed  into  the  oxide  layer  a  space  charge  layer 
develops  across  the  oxide  which  slows  down  the  rate  of  electron  transport.  This 


translates  into  slower  growth  rates  for  the  oxide  film,  after  initial  oxidation.  The 
governing  equation  is  shown  below. 


d  =  k0  In 


(2) 


where  d  is  oxide  thickness,  t  is  time,  t  is  a  time  constant,  and  ko  is  an  initial  oxide 
thickness  that  depends  on  t. 

Experimental  Procedure 

Sample  surface  preparation  is  a  very  important  step  in  oxidation  experiments. 
Surface  roughness,  purity,  and  microstructure  all  impact  the  oxidation  behavior  of  the 
material.  For  repeatable  oxide  thickness  measurements  it  is  necessary  to  create  a  sample 
with  an  ultra  smooth  surface  condition.  Samples  with  a  rough  surface  will  have  greater 
oxidation  because  of  the  increased  surface  area. 

Creation  of  a  smooth  surface  finish  for  indium  initially  proved  to  be  quite 
difficult.  Traditional  mechanical  polishing  of  indium  produces  a  poor  surface  finish 
because  of  indium's  ductility  at  room  temperature.  Silicon  carbide  as  well  as  silica  slurry 
particles  embed  on  indium's  surface.  This  obviously  would  impact  the  oxidation 
behavior,  kinetics,  mechanical  response,  etc. 

To  solve  this  problem  we  adopted  an  electrochemical  polishing  (ECP)  approach. 
To  implement  this  approach  we  built  an  electrochemical  polisher  (seen  in  Fig  1-5  below) 
which  consisted  of  a  rectifier,  anode,  cathode,  and  a  beaker  for  the  electrolyte  solution. 
The  electrolyte  solution  was  a  3:1  ethanol  to  nitric  acid  bath  which  had  to  be  maintained 
at  0°C  because  the  solution  becomes  unstable  at  room  temperature. 

ECP  has  many  benefits  over  traditional  mechanical  polishing.  First,  ECP  only 
takes  a  fraction  of  the  time  mechanical  polishing  takes.  To  mechanically  polish  indium 
to  a  scratch  free  surface  takes  well  over  an  hour  because  of  the  many  steps  involved, 
whereas  ECP  only  takes  2-5  minutes.  Secondly,  ECP  produces  a  much  smoother  finish 
compared  to  mechanical  polishing  because  material  is  removed  from  the  peaks  at  a  much 
faster  rate  than  the  valleys.  Third,  ECP  is  a  stress  free  polishing  technique.  Unlike 
mechanical  polishing  where  residual  stresses  develop  on  the  surface,  ECP  does  no 
mechanical  work  on  the  surface.  Finally,  ECP  produces  a  hygienically  clean  debris  free 
surface  finish.  As  mentioned  earlier,  particles  embedded  in  the  sample  surface  due  to 
mechanical  polishing  will  affect  the  oxidation  behavior. 


(d)  (e)  (f) 

Fig.  1-6.  Optical  and  SEM  images  of  indium  mechanically  polished  (a)  (b)  (c) ,  and 
electrochemically  polished  (d)  (e)  (f). 


The  images  shown  above  in  Fig.  1-6  show  a  side  by  side  comparison  of  indium 
mechanically  polished  and  electrochemically  polished.  As  we  can  see  in  images  (a),  (b), 
and  (c)  mechanical  polishing  produces  a  poor  surface  finish.  Silicon  carbide  particles 
cover  the  sample  surface,  making  repeatable  thickness  measurements  near  impossible. 
However  in  images  (d),  (e),  (f)  we  see  an  ultra  smooth  surface  finish  where  grain 
boundaries  are  exposed. 


Sample  Preparation  and  Temperature  Measurements 


First  the  samples  were  cut  by  razorblade  into  50mm  x  5mm  strips.  Then  these 
strips  were  ECP  in  a  3:1  ethanol  to  nitric  acid  solution  for  5  minutes  at  approximately  4 
volts.  Immediately  after  polishing,  the  samples  were  rinsed  in  distilled  water  for  2 
minutes  and  then  dried  and  stored  in  a  nitrogen  environment.  Before  oxidation  the 
samples  were  cut  into  5mm  x  5mm  squares.  Oxidation  was  carried  out  by  placing  the 
samples  on  a  glass  slide  which  was  placed  on  a  hot  plate.  The  temperature  of  the  solder 
was  monitored  by  a  K-type  thermocouple  which  was  placed  about  1mm  from  the  sample 
surface. 

Thickness  Measurements  through  an  Ellipsometer  (<  40nm) 

After  the  samples  were  thermally  oxidized  thickness  measurements  were  carried 
out  using  NanoFilm’s  EP3  Spectroscopic  Ellipsometer.  Ellipsometry  is  an  especially 
useful  technique  because  of  its  nondestructive  nature  with  little  or  no  special  preparation 
needed  other  than  a  relatively  smooth  sample  surface.  As  demonstrated  in  Fig.  1-7, 
ellipsometry  measures  the  thickness  of  thin  films  by  reflecting  known  polarized  light  onto 
the  sample.  This  polarized  light  reflects  off  both  the  thin  film  and  the  thin  film-substrate 
interface  resulting  in  a  different  polarization  of  the  light.  This  change  in  polarization  is 
analyzed  in  terms  of  s  and  p  light  which  can  be  converted  to  the  familiar  ellipsometry 
variables  ¥  and  A  via  the  following  equation: 

,‘A 

=  tan  y/e  (3) 

,  where  Rp  is  the  reflected  p-light,  and  Rs  is  the  reflected  s-light. 

A  sample  plot  of  A  verses  wavelength  is  shown  in  Fig.  1-8.  In  the  plot  A 
measured  via  ellipsometry  is  compared  to  A  calculated  theoretically  from  the  index  of 
refraction  and  extinction  coefficient  for  indium-indium  oxide  system.  The  mean  square 
error  in  this  plot  is  then  plotted  verses  film  thickness  to  derive  the  film  thickness. 
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Fig.  1-7  Illustration  of  the  ellipsometry  measurement  technique. 
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Fig.  1-8  Plot  of  A  calculated  and  measured  verses  wavelength 


Nanoindentation  (thickness  measurements  >  40nm) 

Nanoindentation  like  other  hardness  tests  involves  the  penetration  of  the  sample 
with  an  indenter.  However,  in  nanoindentation  a  new  depth  sensing  parameter,  contact 
depth  (hc),  is  introduced  (Fig.  9).  The  contact  depth,  the  depth  to  which  the  sample  is  in 
contact  with  the  indenter,  is  measured  by  extending  the  slope  of  the  upper  portion  of  the 
unloading  curve  to  the  x-intercept  as  below: 
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where  e  is  the  geometric  constant. 


(4) 


Fig.  1-9.  Load  (P)  vs.  displacement  (h) 
plot  for  a  typical  nanoindentation 
loading  and  unloading  cycle. 


The  area  function  then  uses  this  contact  depth  to  estimate  the  projected  contact 
area.  Contact  area  can  be  used  to  determine  the  hardness,  H  and  reduced  modulus,  Er  of 
a  sample  from  the  following  equations: 


* 


(5) 


where  P  is  load  and  A  is  the  projected  contact  area,  and: 
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where  S  is  contact  stiffness  and  (3  is  an  indenter  constant.  To  obtain  Young's  modulus 
from  the  reduced  modulus  the  follow  equation  can  be  used: 


1  _  1-t;2  1  -of 

Er~  E  +  E, 

where  subscript  i  indicates  for  diamond  indenter,  and  v  is  Poisson's  ratio. 


(7) 


Technical  Results  and  Data 


Oxidation  and  Reduction 


The  first  oxidation  experiment  was  the  measurement  of  oxide  thickness  verses 
temperature  in  both  air  and  in  an  inert  environment.  All  samples  for  this  experiment  were 
heat  treated  for  two  hours.  A  plot  of  the  results  for  the  samples  oxidized  in  air  can  be  in 
Fig.  1-10. 

Oxidation  of  Pure  Indium  in  Air 
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Fig.  1-10  Thickness  verses  temperature  in  air. 


Several  important  points  can  be  derived  from  this  plot.  First,  there  is  little  or  no 
oxidation  before  the  melting  point  (I56°C).  This  means  indium  can  be  stored  at  elevated 
temperatures  (below  Tm)  without  significant  oxidation.  However,  at  and  above  Tm  there 
is  a  significant  increase  in  oxide  layer  growth  which  shows  asymptotic  behavior  after 
200°C.  So,  from  the  oxidation  point  of  view  there  is  no  significant  difference  if  we 
choose  a  reflow  temperature  of  1 80°C  or  200°C  or  even  220°C.  They  all  show  similar 
oxidation  behavior. 

Conversely,  oxidation  in  an  inert  environment  depends  heavily  on  temperature  in 
that  range.  The  samples  in  the  plot  below  (Fig.  1-11)  were  heat  treated  in  a  glove  box 
with  C>2<  .lppm,  H2O  =  0.1-0.3ppm,  H2=  0.6%.  Here  we  can  see  how  the  oxidation 
behavior  follows  thermodynamic  predictions.  In  Fig.  1-2  we  see  that  at  0.1-0.3ppm 
moisture  content  only  0. 5-0.6%  hydrogen  is  needed  to  create  a  reducing  environment. 
The  chart  below  seems  to  agree  with  this  prediction.  Indium  is  oxidizing  in  the 
temperature  range  from  20°C-160°C.  However  after  160°C  (i.e.  after  crossing  into  the 
reduction  zone)  the  oxide  thickness  seems  to  be  decreasing.  Samples  were  also  heat 
treated  at  200°  C  and  300°  C  at  this  condition  but  measurements  of  these  films  were 
difficult  because  their  shape  became  spherical. 


Inert  Heat  Treatment 
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Fig.  1-11  Thickness  verses  temperature  in  glove  box:  02<  .lppm,  H2O  =  0.1-0.3ppm, 

H2=  0.6% 


We  have  tested  several  different  conditions  (oxidizing,  reducing,  slightly 
oxidizing)  to  confirm  the  validity  of  thermodynamic  oxidation/reduction  map.  These 
data  are  marked  as  symbols  on  the  map,  and  showed  good  correlation  between 
experimental  observations  and  theoretical  predictions  (Fig.  1-12). 
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Fig.  1-12.  Oxidation/reduction  map:  Effect  of  hydrogen  on  1^03  stability  and  transition 
temperature.  Symbols  indicated  experimental  observation  of  indium  samples  at  specific 

conditions. 

Kinetics  of  Oxidation 

The  growth  rate  kinetics  of  indium  was  investigated  in  both  air  and  inert 
environments.  Samples  were  prepared  and  oxidized  from  5  to  120  minutes  at  three 
different  temperatures  in  air  (Fig.  1-13).  After  curve  fitting  all  three  kinetics  curves  to 
the  three  kinetics  mechanisms  described  above,  the  logarithmic  rate  law  showed  the  best 
fit,  all  with  correlation  coefficients  greater  than  0.9.  The  primary  reason  for  such  a  fit  is 
the  asymptotic  behavior  of  the  curves  after  60  minutes.  This  can  be  explained  by  the 
space  charge  layer  limiting  the  flow  of  electrons  after  60  minutes. 


Time  (min) 

Fig.  1-13.  Growth  kinetics  of  indium  in  air  at  145°C,  180°C,  and  220°C. 


Another  supporting  piece  of  evidence  for  logarithmic  oxidation  is  the  initial 
activation  energy  barrier  when  oxygen  is  chemisorbed  onto  the  surface  of  indium. 
Theoretical  calculation  from  the  Rideal  and  Wansbrough-Jones  relation: 

AE  =  0-K  (8) 

where  AE  is  the  activation  energy,  <J>  is  the  work  function,  and  K  is  the  electron  affinity 
of  the  oxygen  atom.  Theoretical  calculation  of  the  activation  energy  for  initial  oxidation 
is  ,52eV.  As  seen  in  the  graph  (Fig.  1-14)  the  experimental  value  of  this  energy  is  ,65eV 
(62.6  kJ/mol),  which  is  quite  consistent  with  the  calculation. 

100 

"c" 

'E 
E 

c 
a 
to 

a: 

<=  10 
o 

to 
■g 
x 

O 


1 

0.00024  0.00025  0.00026  0.00027  0.00028  0.00029 

1/RT  (mol/J) 

Fig.  1-14  Initial  activation  energy  experimentally  determined  to  be  62.6  kJ/mol  (.65eV). 

Oxidation  kinetics  experiments  were  also  carried  out  in  the  glove  box  (Fig.  1-15). 
The  preliminary  results  from  this  experiment  are  shown  in  the  plots  below.  It  seems  there 
is  no  discernable  pattern  in  these  results  since  all  the  values  fall  within  the  range  of  the 
native  oxide  layer.  So  it  would  be  safe  to  assume  the  growth  kinetics  are  very  slow  in 
such  an  environment  (02<  .lppm,  H2O  =  2.3ppm,  H2=  4%). 

Oxidation  Kinetics  at  145°C  (2.3ppm  H20) 
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Fig.  1-15.  Growth  kinetics  of  indium  in  glovebox  (C>2<  .lppm,  H2O  =  2.3ppm,  H2=  4%) 

at  145°C. 


Reduction 


Reduction  of  indium  oxide  from  a  pure  indium  substrate  proved  to  be  quite  a  difficult 
task  due  to  indium  oxides  stability  at  high  temperatures.  We  first  attempted  to  reduce 
indium  oxide  using  argon  plasma  cleaner  to  bombard  the  oxide  with  ions.  Preliminary 
results  from  this  experiment  shown  in  Fig.  1-16,  where  an  etching  rate  of  5.3  nm/hr  was 
observed  at  270  watts. 
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Fig.  1-16.  Reduction  of  indium  oxide  via  an  argon  plasma  cleaner. 

Reduction  was  also  attempted  by  adding  hydrogen  to  the  glove  box  environment 
to  make  indium  oxide  unstable  relative  to  indium  at  oxygen  partial  pressures  much 
greater  than  the  equilibrium  oxygen  partial  pressure.  Figure  1-17  shows  a 
thermodynamic  map  where  samples  are  oxidized  at  conditions  left  of  the  equilibrium 
curve,  and  reduced  at  conditions  to  the  right  of  the  curve.  Many  experiments  were 
preformed  in  the  reduction  zone  with  no  noticeable  reduction  of  the  oxide  layer.  Only  a 
harsh  reducing  environment  showed  a  reduction  in  the  oxide  layer.  With  the  temperature 
set  to  350°C,  H2  =  4.5%,  O2  <  .Ippm,  and  H2O  =  7ppm  the  oxide  layer  was  reduced  from 
21nm  to  15nm  after  6  hours,  yielding  an  etching  rate  of  lnm  per  hour,  as  shown  in  Fig. 

1 8.  This  clearly  demonstrates  the  need  for  very  low  moisture  levels  so  the  temperature 
can  be  reduced  into  an  acceptable  range. 


Reduction  Map  At  7ppm  Moisture 
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Fig.  1-17  Thermodynamic  map  of  indium  oxide  reduction  at  7ppm  moisture. 


Hydrogen  reduction  of  Indium  Oxide 


Fig.  1-18.  Reduction  of  indium  oxide  via  a  4.5%  hydrogen  gas. 


Figure  1-19  shows  the  sample  coupons  before  reduction  experiments  (i.e.,  21  nm 
thick  oxide  on  indium)  and  after  reduction  experiments  (i.e.,  15  nm  thick).  After  being 
exposed  to  reducing  environment,  the  coupon  became  more  rounded,  indicating  that 
oxide  layer  gets  somewhat  removed. 


Fig.  1-19  (a)  Sample  before  reduction  (b)  sample  after  reduction. 


Thickness  measurements  via  Nanoindentation 


Dynamic  nanoindentation  superimposes  a  sinusoidal  load  on  top  of  a  quasi-static 
load  to  obtain  contact  stiffness  measurements  continuously  as  a  function  of  displacement. 
Given  the  material  system  indium-oxide/indium  there  is  an  obvious  change  in  the  contact 
stiffness  as  the  indenter  travels  past  the  film  into  the  substrate.  This  allowed  us  to 
estimate  the  thickness  of  several  films  grown  at  higher  temperatures.  Figure  1-16  shows 
both  the  dynamic  and  quasi-static  nanoindentation  plots  where  there  is  a  sudden 
displacement  burst  as  the  indenter  travels  past  the  oxide.  Figure  1-20  (a)  and  (b)  show 
this  behavior  at  a  displacement  around  100  nm  whereas  Fig.  1-20  (c)  and  (d)  show  this 
behavior  at  around  800  nm. 


Indium  Melted  at  500°C  For  Ihour 
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Indium  Sample  Melted  at  600'C 
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Fig.  1-20  Dynamic  and  quasi-static  nanoindentation  results  for  indium  heat  treated  at 

500°C  and  600°C  respectively. 

Materials  Characterization 


NanoDMA  Estimation  of  Elastic  Modulus  on  Indium 


NanoDMA  Estimation  of  Hardness  for  Indium 
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Indium  Melted  at  1000°C 


Indium  Melted  at  1000°C 


a  ••  •  a  a**  ••  •  «.  »**•  ^ 

%xx  X  X*  ***  X**  ***/*«T x  If. 

^  ^  t  -  *■  a-  » r  it  f v : .  *  ^  ^ 


100  WO  140  160 


*°  700  ?4°  w>  780  150  170  190  210  230  250  270 

Displacement  (nm)  _.  >  .  .  . 

r  *  Displacement  (nm) 

(b)  (d) 


Fig.  1-21  (a)  elastic  modulus  of  indium;  (b)  hardness  of  indium;  (c)  elastic  modulus  of 
indium  oxide;  (d)  hardness  of  indium  oxide,  (note  that  thick  indium  oxides  were  grown 

from  melting  at  high  temperatures) 


Materials  characterization  of  both  indium  and  indium-oxide  where  carried  out 
using  nanoindentation,  optical  microscopy,  AFM,  SEM,  and  XRD.  Results  from  the 
mechanical  characterization  using  dynamic  nanoindentation  are  shown  below. 

The  elastic  modulus  of  indium  was  measured  to  be  12.4  ±.34  GPa,  and  the 
hardness  was  39.7  ±  1.2  GPa.  The  modulus  of  indium  oxide  was  measured  to  be  105.85 
±  7.6  GPa  and  the  hardness  6.64  ±  .96  GPa  (Fig.  1-21).  Microstructural  characterization 
of  indium  oxide  annealed  at  1000  °C  was  preformed  using  AFM.  From  Fig.  1-22  we  can 
see  that  the  microstructure  is  composed  of  oval  grains  200-300  nm  in  size.  The  grains  of 
indium  seen  earlier  are  much  larger  2-4  mm. 


Fig.  1-22  AFM  images  of  indium  oxide  annealed  at  1000  °C. 

X-ray  Diffraction  (XRD)  Analysis 

The  crystal  structure  of  the  TGO  grown  on  indium  samples  was  analyzed  using 
X-ray  diffraction.  Samples  below  and  above  the  melting  temperature  were  tested  to 
determine  any  changes  in  the  crystal  structure  of  indium  oxide.  Below  are  the  resulting 
X-ray  intensity  verses  20  plots.  It  can  be  seen  from  these  plots,  that  below  the  melting 
temperature  no  crystalline  phase  of  indium  oxide  was  detected. 


Fig.  1-23  XRD  plot  of  intensity  verses  2theata  for  the  sample  annealed  at  153°C  for  4 
hours.  The  peaks  labeled  are  from  pure  indium  metal  (substrate). 


Figure  1-23  represents  the  XRD  plot  of  indium  annealed  at  153°C  for  4  hours. 
Here  the  peaks  only  correspond  to  the  crystal  structure  of  the  pure  indium.  Therefore,  it 
is  believed  that  indium  oxide  grown  below  the  melting  temperature  is  of  an  amorphous 
nature.  Conversely,  when  indium  was  annealed  above  the  melting  temperature  a  cubic 
indium  oxide  crystal  structure  was  detected.  Below  in  Fig.  1-24  the  peaks  correspond  to 
a  cubic  crystal  structure  for  indium  oxide  when  the  sample  was  annealed  for  four  hours  at 
200°C.  When  treated  above  the  melting  temp  the  sample  begins  to  flow  resulting  in  a 
rough,  deformed  surface.  This  may  also  explain  the  excess  “noise”  in  the  plot. 


Fig.  1-24  XRD  plot  of  intensity  verses  2-theata  for  the  sample  annealed  at  200C  for  4 
hours.  The  peaks  labeled  are  from  indium  oxide  (film). 


Conclusion 

There  are  several  conclusions  that  can  be  made  regarding  the  creation  of  a  solder 
reflow  processing  window.  From  the  thermodynamics  point  of  view  the  four  main 
parameters,  temperature,  oxygen  content,  moisture  content  and  hydrogen  content  can  be 
precisely  controlled  to  produce  a  reducing  or  inert  ( reducing  but  very  slow  etching  rate) 
environment.  Our  calculations  revealed  that  moisture  is  a  critical  parameter  in  the 
creation  of  a  reducing  environment.  Reducing  the  moisture  to  the  lowest  possible  levels 
shifts  the  equilibrium  curves  to  the  left  which  lowers  the  temperature  and  hydrogen 
content  needed  to  create  the  reducing  environment.  We  also  observed  that  oxygen 
content  is  not  a  critical  parameter  for  producing  a  reducing  environment.  Shifting  the 
oxygen  content  from  lppm  to  500ppm  will  not  shift  the  location  of  thermodynamic 
equilibrium. 

With  regards  to  the  kinetics  investigation,  we  found  that  the  growth  rate  is 
controlled  by  a  logarithmic  relationship.  In  the  growth  rate  equation  the  parameters  t  and 
ko  reveal  the  initial  onset  of  oxide  growth  ko  in  timex.  If  the  reflow  time  can  be  held 
below  x  significant  oxidation  can  be  avoided.  We  also  saw  that  that  there  is  little  or  no 
oxidation  samples  when  samples  were  heat  treated  in  the  glove  box  below  the  melting 


temperature.  This  means  samples  can  be  stored  in  the  glove  box  with  no  significant 
additional  oxidation  (other  than  native  oxide). 

As  for  the  reduction  of  indium  oxide  in  a  hydrogen  environment,  we  found  the 
kinetics  to  be  very  slow  in  regions  close  to  thermodynamic  equilibrium.  Only  very  harsh 
conditions  (T  >  350°C)  reduced  the  oxide.  However  by  further  reducing  the  moisture 
content,  the  service  temperature  should  also  be  reduced  into  an  acceptable  range.  In  this 
case,  pre-reflow  conditions  for  an  extended  period  prior  to  joining  may  enable  to  remove 
initial  oxides  formed  on  the  surface  of  indium.  More  experiments  are  warranted  to 
quantitatively  verify  the  reduction  kinetics.  Using  the  attached  thermodynamic  program 
(see  Appendix  A)  will  guide  the  search  for  proper  reflow  environments. 

Appendix  A. 

Thermodynamic  Reflow  Code  (MATLAB  Program): 

%Harry  Schoeller 
%Thermodynamic  Program 

%This  program  determines  whether  your  environment  will  oxidize  or  reduce 
%indium  solder 

Y  =  input('Welcome  To  The  Indium  Solder  Reflow  Environment  Calculator:  Run  (1)  or 
Exit  (2):  '); 
while  Y  ==  1, 

Z  =  inputC\n Which  parameter  would  you  like  to  hold  constant  pH20  (1),  pH2  (2),  pH20 
and  pH2  (3):  ’); 
ifZ  ==  1, 

P  =  input('Enter  the  total  pressure  of  the  system  (atm):  '); 
pH201=  input('Enter  the  pH20  (ppm):  ’); 
p021=  input('Enter  the  p02  (ppm):  '); 

T1  =  input('Enter  the  reflow  temperature  (C):  ’); 

p02 1  =  p02 1  *  (10A-6); 

pH20  =  pH20 1  *(  1 0A-6); 

T  =  T1  +273; 
if  T  <=  430, 

G  =  183300  -  151.5*T; 
else 

G  =  189300-  168.1  *T; 
end 

H  =  (pH20/(exp(G/(-8.3 1 44*T))A(  1/3)))*  1 00; 
if  H  >  100, 

fprintf(’lndium  oxide  cannot  be  reduced  with  this  moisture/temperature 
combination,  please  increase  temperature  or  reduce  moisture  content’); 
else 

fprintf('Your  environment  needs  at  least  %4.4f  percent  hydrogen  gas  to  reduce 
indium  oxide',  H); 
end 


Y  =  input('\n Welcome  To  The  Solder  Reflow  Environment  Calculator:  Run  (1)  or 
Exit  (2):  ’); 

elseifZ  ==  2 

P  =  input('Enter  the  total  pressure  of  the  system  (atm):  ’); 
p021=  input('Enter  the  p02  (ppm):  '); 

H=  input('Enter  the  percentage  of  hydrogen  gas:  ’); 

T1  =  input(’Enter  the  reflow  temperature  (C):  '); 

T  =  T1  +  273; 
if  T  <=  430, 

G  =  183300  -  151.5*T; 
else 

G=  189300-  168.1  *T; 
end 

pH20  =  ((H/100)*  (exp(G/(-8.3 1 44*T))A(  1  /3))); 

p021  =  p021  *  (10A-6); 

ppm  =  pH20/(.00001); 

fprintf('Your  environment  should  contain  less  than  %4.4f  ppm  H20  to  reduce  indium 
oxide’,  ppm); 

Y  =  input('\nWelcome  To  The  Solder  Reflow  Environment  Calculator:  Run  (1)  or 
Exit  (2):  '); 

elseifZ  ==  3 

P  =  input('Enter  the  total  pressure  of  the  system  (atm):  ’); 
p021=  input('Enter  the  p02  (ppm):  '); 
pH201=  input('Enter  the  pH20  (ppm):  '); 

H=  input('Enter  the  percentage  of  hydrogen  gas:  ’); 
pH20  =  pE120 1  *(  1 0A-6); 

p021  =  p021*  (10A-6); 

HI  =  H/100; 

T  =  (-22767.7/(log((pH20/H  1  )A3)-20.2)); 

T2  =  T  -  273; 

fprintf('Your  reflow  temperature  should  be  at  least  %4.2f  C  to  reduce  indium  oxide', 
T2); 

Y  =  input('\n Welcome  To  The  Solder  Reflow  Environment  Calculator:  Run  (1)  or  Exit 
(2):  '); 
end 
end 


1-2.  Wettability/Solderability 


Sample  preparation 

Indium  solderability  was  assessed  by  wetting  angle  measurement  and  joint  strength 
measurement  with  a  single  lap  shear  test.  Indium  oxide  thickness  was  taken  of  an 
experimental  parameter  to  affect  the  solderability.  The  indium  oxide  was  grown  at 
different  temperatures(25  °C,  145  °C,  160  °C  and  200  °C)  for  2  hours  on  a  hot  plate.  The 
measured  oxide  thickness  is  listed  in  the  table  1-1. 


Table  1-1  Indium  oxide  thickness  vs.  tern 

aerature  (2  hours  annealing) 

V^p(°C) 

Oxiarv. 

thicknesses^ 

25 

145 

160 

200 

AVG  (nm) 

4.0 

5.2 

13.9 

16.9 

Before  the  indium  solderability  test,  the  shape  of  indium-bonding  substrate(6.5  by  28 
mm)  and  size  of  bonding  pad  were  manufactured  on  the  silicon  wafer  by  micro¬ 
fabrication  process  in  Figure  1-25.  Au  and  titanium  were  then  deposited  on  the  several 
silicon  wafers.  In  the  region  except  for  indium  bond  area(3.5  by  3.5  mm),  Au  and 
titanium  were  clearly  removed  by  chemically  dissolving  method  in  Figure  1-26.  The 
deposited  thickness  of  Au  layer  is  0.5  urn,  and  that  of  titanium  is  0.04  urn.  Titanium  was 
used  to  improve  the  adhesion  of  indium  to  the  silicon  substrate,  and  Au  was  used  to 
protect  the  titanium  from  oxidation  in  air. 


Au  layer 

Tilanium  layer 


Fig.  1-25  Silicon  substrate  design  Fig.  1-26  Detailed  size  of  silicon  substrate 


Wetting  angle  measurement 

The  indium  for  wetting  angle  measurements  was  melt  on  the  Au-deposited  silicon 
wafer  using  hot  plate(200  °C  for  2  min.)  in  an  air  and  in  a  glove  box.  The  atmosphere 
condition  in  a  glove  box  was  controlled  by  the  concentration  of  moisture,  nitrogen  and 
hydrogen  gas.  The  glove  box  can  be  operated  at  H2O  <  1  ppm,  O2  <  0.1  ppm,  and  H2  <  5 


%.  After  the  wetting  was  formed,  wetting  angles  were  measured  by  Wyco  surface  profiler 
and  were  determined  at  the  interface  of  indium  and  the  Au  layer  on  the  substrate. 

Pure  indium(with  4  nm  oxide  thickness)  wetting  shapes  at  an  ambient  condition  and 
in  a  controlled  environment(C>2<  0.1  ppm,  H20<  1  ppm,  and  H2=l  .5  %)  are  compared 
with  each  other  in  Figure  1 -27(a)  and  (b).  The  wetted  indium  in  an  air  ambient  condition 
is  just  melted  on  the  Au  layer  in  Figure  1 -27(a),  but  the  wetted  indium  in  a  controlled 
glove  box  environment  shows  the  pronounced  spreading  after  melting  on  a  hot  plate  at 
200  °C,  as  shown  in  Figure  1 -27(b).  After  all,  the  indium  wetted  in  a  controlled  glove  box 
environmental  condition  shows  better  wetting  than  that  in  an  air  ambient  condition. 


Fig.  1-27  Wetting  shapes  of  indium 

Figure  1-28  shows  wetting  angle  measurements  with  initial  indium  oxide  thickness 
quantitatively.  In  a  case  of  4  nm  thickness  indium  oxide,  indium  wetting  angle  on  Au 
layer  is  about  1 5  degree.  The  wetting  angle  was  increased  to  30  degree  at  1 7  nm 
thickness  indium  oxide. 


Initial  indium  oxide  thickness,  nm 


Fig.l-28  Wetting  angle  variation  with  initial  oxide  thickness  in  a  controlled  environment 


Joint  strength  measurement 


Indium  was  reflowed  on  the  Au  layer  between  a  pair  of  silicon  substrates  for  single  lap 
shear  joint  test  at  different  environments.  The  test  geometry  is  shown  in  Figure  1-29.  First 
of  all,  indium  reflow  for  bonding  was  performed  on  the  hot  plate  using  a  specially 
designed  fixture(Figure  1-30)  in  an  air.  The  reflow  condition  of  indium  is  around  230  °C 
for  10  minute  in  an  air  ambient  condition. 


indium 


Unit:  mm 


Fig.l-29  Test  geometry  for  joint  shear  test 
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Fig.1-30  Design  of  test  fixture  for  joint  shear  test 

Figure  1-31  shows  a  real  joint  test  sample  with  4  nm  oxide  thickness  before  joint 
test.  Indium  was  embedded  and  reflowed  to  each  Au  layer  between  two  silicon  substrates. 
After  performing  indium  joint  failure  test  by  single  lap  shear,  the  failure  was  dominated 
by  interfacial  mode  failure  between  indium  and  the  substrate,  as  shown  in  Fig.  1-32.  The 
40  time  magnified  view  of  Fig. 1-32  shows  a  lot  of  void  traces  formed  at  the  interface 
between  indium  and  the  substrate  in  Figure  1-33. 


Fig.1-31  Real  bonded  sample  for  joint  test 


Fig.1-32  Interfacial  failure  between  indium  and  bond  substrate 


Fig.1-33  Magnified  view(40  X)  of  indium  failure  site  in  Figure  8 


Figure  1  -34  shows  one  of  representative  curves  which  were  tested  for  indium  with  a 
4  nm  oxide  thickness  in  an  air  reflow  condition.  The  reflow  temperature  of  indium  is  230 
°C.  The  indium  joint  strength  by  a  single  lab  shear  is  1 .35  Mpa  with  standard  deviation  of 
0.2  Mpa  at  a  strain  rate  of  0.5  mm/sec. 

Figure  1-35  shows  the  joint  strength  variation  with  the  change  of  initial  oxide 
thickness.  As  initial  oxide  thickness  is  increased,  indium  joint  strength  is  decreased  to 
below  1 .0  Mpa  due  to  the  existence  of  void  defect  and  of  much  thicker  indium  oxide. 

Environmental  condition  in  a  glove  box  can  be  controlled  by  the  concentrations  of 
hydrogen  and  moisture  in  Figure  1-36,  which  is  drawn  by  thermodynamic  calculation  for 
indium  oxidation  phenomena.  An  oxidizing  condition  is  reached  below  each  line 
corresponding  to  the  moisture  concentration.  On  the  above  line,  the  condition  will  be  in  a 
reduced  environment.  On  the  line,  the  atmosphere  is  in  an  inert  environment. 


Fig.1-34  Indium  joint  strength  with  4  nm  initial  oxide  thickness  in  an  air  reflow 
condition(230  °C  for  more  than  5  minutes) 


Fig.1-35  Indium  joint  strength  with  initial  oxide  thickness  in  an  air  reflow  condition  (230 

°C  for  more  than  5  minutes) 
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Fig.l-36  Moisture  effect  on  the  environmental  condition  in  a  glove  box 

Figure  1-37  shows  indium  joint  strengths  with  initial  different  oxide  thicknesses.  The 
joint  strength  was  dropped  with  the  increase  of  initial  oxide  thickness.  Those  tested 
samples  were  made  by  using  specially  designed  fixture(Fig.3.6)  at  controlled  oxidizing 
glove  box  conditional)  <  6  ppm,  02  <  0.1  ppm,  and  H2  <  4  %).  The  reflow  condition  for 
indium  bonding  on  the  substrate  was  160  ~  170  °C  for  5  minutes.  The  joint  strength  is 
measured  by  tensile  loading  at  2  mm/sec  strain  rate. 


Displacement,  mm 

Fig.l-37  The  comparison  of  indium  joint  strengths  with  initial  different  oxide  thickness 
at  controlled  environmental  glove  box  condition 

Figure  1-38  is  plotted  to  compare  indium  joint  strength  made  at  a  controlled 
environment  with  that  of  an  air  ambient  environment.  There  is  a  big  difference  of  joint 
strength  depending  on  the  environmental  conditions  that  the  samples  were  prepared. 
Those  testing  results  are  obtained  by  the  samples  which  are  made  with  the  same  specially 
designed  fixture.  The  fixture  was  designed  to  apply  hot  pressure  on  indium  when  it  was 


melt.  In  order  to  obtain  the  indium  joint  strength  on  Au  layer  without  pressure,  indium 
was  put  to  be  melt  on  the  deposited  area  of  Au  on  the  silicon  substrate  without  any 
pressure,  and  then  tested  by  die  shear  load  at  2  mm/sec  strain  rate  using  Dage  bond  tester. 
The  die  shear  load  results  are  also  added  to  the  Fig.  1-38. 
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Fig.l-38  The  variation  of  indium  joint  strength  at  different  environmental  condition  and 

indium  bond  condition 

Figure  1-39  shows  the  trace  of  indium  bond  site  after  joint  failure  test.  The  Au  and 
titanium  were  observed  at  one  side  of  bond  substrate  and  the  indium  was  remained  at  the 
other  side  substrate. 


Fig.1-39  View  of  bonded  surface  after  joint  failure  at  the  controlled  environmental  glove 

box  condition 


Sub  Task  2:  Thermo-mechanical  Analyses 


Performer  -  Bahgat  Sammakia  /  Bill  Infantolino  /  Harita  Machiraju 
Organization  -  JEEC  -  Binghamton  University 

Introduction 

The  objective  of  the  thermal  analysis  was  to  obtain  the  temperature  distribution  of  the 
complete  S&A  chip  package  which  consists  of  the  SOI  (silicon  on  insulator)  chip  bonded 
to  the  kovar  carrier.  The  boundary  conditions  were  estimated  based  on  the  models 
provided  by  Indian  Head.  The  model  prepared  based  on  the  design  information  received, 
was  called  the  base  case  model.  The  temperature  distribution  for  the  base  case  model  was 
obtained.  A  set  of  parametric  studies  and  design  changes  were  studied  in  order  to  reduce 
the  power  requirements  for  the  actuators.  The  base  case  model  was  used  as  a  reference 
for  all  the  parametric  studies  performed.  A  computational  fluid  dynamics  model  was  also 
prepared  to  determine  the  pressure  rise  due  to  the  air  temperature  increase  within  the 
sealed  package. 

The  S&A  chip  consists  of  three  sets  of  micro  actuators  -  arming,  flow  sensor  and  g- 
sensor  actuators,  on  its  device  layer  (as  shown  in  Figure  2-1).  The  micro  actuators,  when 
powered,  provide  the  necessary  displacement  to  perform  a  specific  function,  e.g.,  move 
an  optical  fiber.  The  micro-actuator  is  essentially  a  beam  which  deflects  when  heated  due 
to  thermal  expansion. 

A  three  dimensional  model  was  prepared  to  perform  the  thermal  analysis.  The  initial 
model  was  a  single  beam  model  which  was  done  to  gain  an  understanding  of  the 
temperature  distribution  in  the  beam  and  the  effects  of  various  types  of  heat  transfer  on 
the  beam  temperature.  It  was  also  done  to  compare  the  numerical  results  to  known 
analytical  solutions  with  similar  boundary  conditions. 
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Fig.  2-1:  S&A  chip  with  the  three  sets  of  actuators 


2-1.  Thermal  analysis  of  a  micro  beam 


Description  of  model: 


Adiabatic 

T  =  300K 

condition 

\ 


Conduction  through  stagnant  air  gap 


Initial  Condition  =  300K 


Fig.  2-2:  One  half  of  the  micro-cantilever  beam  with  symmetry  loading  conditions 


The  beam  in  the  initial  model,  shown  in  Figure  2-2,  was  considered  to  be  attached 
to  one  anchor  at  each  end.  Since  these  anchors  are  relatively  large  masses  compared  to 
the  beam  itself  and  the  time  for  the  analysis  is  relatively  short,  the  anchors  were  assumed 
to  be  at  a  constant  temperature.  The  power  input  to  the  actuator  was  assumed  to  generate 
heat  uniformly  within  the  whole  beam.  The  heat  generation  is  not,  however,  extended  to 
the  anchors.  Only  half  of  the  beam  was  modeled  due  to  symmetry. 

The  length  of  the  beam  modeled  was  5500  microns,  the  height  of  the  beam,  125 
microns  and  the  width  of  the  beam,  10  microns,  based  on  the  inputs  provided  at  that  time. 
The  numerical  software  used  for  performing  the  analysis  was  ANSYS  and  the  element 
type  used  for  the  beam  was  SOLID  70  (3D  Thermal  element).  The  material  properties  for 
this  analysis  are  shown  in  Table  2. 


Material 

Value 

Thermal  Conductivity  (W  m"1  K'1) 

Silicon 

145 

Air 

0.03 

Emissivity 

Silicon 

0.7 

Density  (kg  m'3) 

Silicon 

2330 

Heat  Capacity  (J  kg^K'1) 

Silicon 

729.614 

Table  2-1 :  Material  properties  used  for  tl 

lis  analysis 

Mesh  Convergence  Study: 


The  conduction  through  the  stagnant  air  gap  at  the  base  of  the  beam  results  in  the 
removal  of  a  significant  amount  of  heat.  However  the  temperature  gradient  along  the 
height  of  the  beam  due  to  this  effect  was  small.  The  gradient  was  high  toward  the  ends  of 
the  beam.  Therefore,  the  number  of  elements  in  the  region  close  to  the  constant 
temperature  end  is  higher  than  in  the  center  of  the  beam.  The  results  of  the  mesh 
convergence  study  are  shown  in  Figure  3.  With  the  help  of  this  study,  the  optimum 
number  of  elements  to  be  used  for  the  analysis  was  determined. 


Mesh  Convergence  Study 


600 


Comparison  of  numerical  model  temperature  profile  to  analytical  solution: 

The  first  case  was  modeled  with  boundary  conditions  which  allowed  comparison 
to  a  known  analytical  solution  for  verification  of  the  analysis.  Only  half  of  the  beam  was 
modeled  since  the  beam  and  loading  conditions  were  symmetric.  The  beam  was  subjected 
to  a  fixed  temperature  at  one  end  and  an  adiabatic  condition  to  simulate  the  symmetry  at 


the  other  end.  A  power  input  of  5  Watts  was  applied.  A  transient  analysis  was  performed 
and  the  maximum  temperatures  were  determined  for  various  times. 

The  maximum  temperatures  reached  in  the  system  are  shown  in  Figure  4.  Note 
that  many  temperatures  are  above  the  melting  point  of  silicon  and  are  calculated  from  the 
finite  element  model  only  for  comparison  to  the  analytical  model. 


Maximum  Temperature  Vs  Time 
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Fig.  2-4:  Comparison  of  analytical  and  numerical  solution  for  single  beam  model  -  fixed 

temperature  at  one  end. 

The  results  for  the  above  case  were  verified  using  the  analytical  solution  with 
the  same  loading  conditions.  The  expression  for  the  temperature  distribution  in  a  semi¬ 
infinite  solid  with  a  fixed  temperature  condition  at  one  end  and  an  adiabatic  condition  at 
the  other  end  [1]  is  given  by 


where  a  =  Initial  temperature  (K) 

a  =  Thermal  diffusivity  (sec  m'2) 
t  =  Time  (seconds) 

A0  =  Heat  generation  rate  (W  m'3) 
k  =  Thermal  conductivity  (W  m'1  K"') 
x  =  Position  along  the  beam  (m) 


Figure  2-5  shows  good  agreement  between  the  analytical  solution  and  the 
finite  element  model.  Calculated  temperatures  shown  in  Table  2-2  at  various  positions 
along  the  beam  are  compared  to  the  model  results.  For  example,  when 

a=300K 

a  =k/pC  =  145  /  (2330  *  729.614)  =  8.5294e-5  m2/sec 
t  =  10  milliseconds 
A0  =  5 W/volume  =  1 .0 1 0 1  e  1 2  W/m3 
K=  145  Win'  K'1 


The  results  show  that  the  error  is  approximately  1  %  for  a  range  of  nodal  positions. 


Nodal  Position(microns) 

ANSYS(K) 

Analytical(K) 

Error  % 

100 

697.381 

705.5165 

1.153127 

200 

1338.56 

1354.8 

1.198701 

400 

2430.89 

2462.2 

1.271627 

600 

3355 

1.322206 

800 

|  4016.49 

4071.1 

1.341406 

Table  2-2:  Base  case  results  -  comparison  with  analytical  solution. 


2-2  Thermal  analysis  of  the  S&A  chip  package 
Description  of  the  package/model: 

The  S&A  chip  is  prepared  on  an  SOI  (silicon-on-insulator)  wafer.  This  consists  of  a 
device  layer  and  a  handle  layer  with  an  insulator  in  between.  The  actuators  are  fabricated 
using  DRIE  on  the  device  layer  (125  microns).  The  handle  layer  is  a  400  micron  thick 
layer  of  silicon.  The  two  layers  are  separated  by  a  thin  layer  of  silicon-dioxide  (2 
microns).  The  S&A  chip  consists  of  three  actuators  -  the  flow  sensor,  arming  and  the  g- 
sensor  actuator  on  the  device  layer.  The  device  layer  is  attached  to  a  top  chip  with  indium 
solder  joints  (10  micron  standoff).  The  handle  layer  is  attached  to  the  base  of  the  kovar 
carrier  with  a  gold-tin  eutectic  bond.  A  clear  picture  of  the  actual  package  is  shown  in 
Figure  2-5.  The  carrier  is  attached  to  the  aluminum  enclosure  through  a  circuit  board 
sandwiched  between  thermal  pads.  The  enclosure  is  then  bolted  to  the  torpedo  wall.  The 
placement  of  the  package  in  the  enclosure  was  determined  from  the  models  (Pro-E) 
provided  by  Indian  Head.  The  S&A  chip  package  and  its  surroundings  are  shown  in 
Figure  2-6. 

A  3-D  finite  element  model  was  prepared  for  the  thermal  analysis  of  the  S&A  chip 
package.  The  software  used  for  the  analysis  was  ANSYS  9.0.  The  model  prepared  was 
called  the  base  case  model  and  was  used  as  a  reference  for  several  parametric  studies. 


/  +  Thermal  pad 

Carrier 

Fig.  2-7  :  Section  view  of  the  package 


Based  on  the  inputs  provided,  the  base  case  model  of  the  S&A  chip  package  included  the 
following  parts. 

•  Lid  -  Kovar 

•  Carrier  -  Kovar 

•  Top  Chip- Silicon 

•  S&A  Chip  package  -  Silicon 

•  Solder  joints  (Top  chip)-  lOOIn 

•  Chip  to  carrier  bond  -  AuSn 


•  Thermal  Pad  (Thermaflo  TCP  -  8020TPL220) 

•  Aluminum  fuze  body 

•  Aluminum  fuze  cap 

(Dimensions  taken  from  AUTOCAD  and  Pro-E  models  sent  by  1HD) 

The  section  view  of  the  package  can  be  seen  in  Figure  2-7.  The  2  micron  layer 
gap  between  the  device  layer  and  handle  layer  of  the  SOI  wafer  was  included  in  the 
model  in  the  form  of  a  contact  resistance.  The  effect  of  introducing  contact  elements 
instead  of  physical  elements  was  studied  using  2-D  models.  This  was  mainly  done 
because  of  meshing  restrictions.  The  results  showed  that  the  difference  in  the  temperature 
profile  between  the  two  cases  was  less  than  0.5  %.  The  thermal  conductance  of  air  was 
used  as  the  characteristic  property  for  the  contact  elements  below  the  actuators  and  the 
thermal  conductance  of  silicon-dioxide  was  used  for  the  elements  below  the  rest  of  the 
device  layer.  The  thermal  conductance  of  air  also  included  the  effect  of  conduction  from 
the  sides  of  the  beam  calculated  in  the  form  of  a  shape  factor.  The  other  parts  of  the 
package  modeled  as  contact  elements  were  the  chip  to  carrier  bond  (Au-Sn),  thermal  pad, 
section  of  board  zero  and  the  conduction  of  heat  through  the  air  gap  between  the  S&A 
chip  and  the  top  chip. 


Fig.  2-8:  Device  layer  of  S&A  chip  with  the  push  blocks 


A  push  block  was  included  in  the  center  of  the  actuator  beams.  The  heat 
generation  was  applied  only  in  the  beam  portion  of  the  actuator  and  not  in  the  push 
blocks.  This  is  because  the  push  block  offers  very  low  resistance  to  the  current  passing 
through  it  compared  to  the  rest  of  the  beam  due  to  the  increase  in  cross-sectional  area  in 
that  region.  This  results  in  very  little  joule  heat  generation  in  the  push  block.  This  was 


again  verified  using  a  simple  3-D  model  of  the  actuator  alone  (beams  with  the  push 
block)  and  applying  a  voltage  difference  across  its  ends.  The  joule  heat  generation  in  the 
push  block  was  negligible.  The  push  blocks  in  the  device  layer  are  shown  in  Figure  2-8. 


The  kovar  carrier  is  attached  to  a  circuit  card  with  a  thermal  pad  between  these 
two  components.  The  circuit  card  is  mounted  on  an  aluminum  enclosure  using  the  same 
material.  The  thermal  resistance  for  the  pads  and  the  card  are  included  in  the  model  using 
a  series  of  contact  resistances.  The  aluminum  enclosure  is  bolted  to  the  torpedo  wall 
(seen  in  Figure  2-9)  which  is  subjected  to  a  forced  convection  boundary  condition.  The 
geometry  of  the  aluminum  enclosure  and  torpedo  wall  are  not  modeled  in  detail  but  are 
approximated  in  a  thermally  equivalent  manner.  In  reality,  the  aluminum  fuze  body  is 
made  up  of  two  parts  as  shown  in  Figure  2-10.  There  is  no  interface  material  between  the 
two  parts  of  the  fuze  body.  Due  to  this  reason,  the  thermal  contact  conductance  between 
these  parts  may  be  very  low.  Considering  only  the  lower  portion  of  the  fuze  body  gives 
the  worst  case  condition.  The  fuze  cap  in  the  model  is  not  the  same  size  as  in  the  Pro-E 
models.  Since  the  convection  boundary  condition  is  applied  on  the  outer  surface  of  this 
fuze  cap,  the  larger  area  has  been  simulated  by  multiplying  the  heat  transfer  coefficient 
by  a  factor  equivalent  to  the  ratio  of  the  actual  surface  area  of  the  fuze  cap  to  the  surface 
area  modeled  for  the  analysis.  All  of  these  components  have  been  included  in  the  model 
to  provide  a  representative  boundary  condition.  The  material  properties  used  for  the 
model  are  shown  in  Table  2-3  and  the  dimensions  of  all  the  parts  of  the  model  are  shown 
in  Table  2-4. 


Fig.  2-9:  Complete  package  mounted  on  simulated  Fig.  2-10:  Fuze  cap  and  fuze  body 
aluminum  enclosure  models  (Pro-E  models  -  IHD) 


Thermal  Conductivity 
(Wm-1K-1) 

Specific  heat 
(J  kg-1  K-l) 

Density 
(kg  m-3) 

Silicon 

T(K)  k(Wm-lk-l) 

300  128.20 

320  117.10 

340  108.80 

350  105.90 

400  89.13 

500  67.46 

600  53.73 

700  44.32 

800  37.51 

1000  28.39 

1200  22.61 

729.614 

2330 

Air 

300  0.02624 

350  0.03003 

400  0.03365 

450  0.03707 

500  0.04038 

550  0.0436 

600  0.04659 

650  0.04953 

700  0.0523 

750  0.05509 

800  0.05779 

100  Indium 

83 

243 

7310 

Kovar 

16 

418 

7750 

180 

898.7 

2700 

AuSn  (Bond) 

57 

- 

Glass  (Corning7059) 

1 

- 

- 

Thermaflo  TCP  - 
8020TPL220  pad 

6 

“ 

Table  2-3:  Material  properties  of  the  different  parts  of  the  S&A  chip  package 


Kovar  carrier-inner  walls 

17.7 

17.7 

3.11 

Kovar  carrier  base 

24.4 

2.45 

Gold-Tin  eutectic 

12.35 

12.35 

25.4 

Table  2-i 

:  Dimensions  of  dif 

erent  parts  of  the  package 

Boundary  conditions: 

The  actuators  are  subjected  to  a  uniform  heat  generation.  As  mentioned  before,  the 
heat  generation  is  applied  only  in  the  beam  and  not  in  the  push  block. 

The  torpedo  wall  is  subjected  to  forced  convection.  The  value  of  the  heat  transfer 
coefficient  used  in  the  model  is  1000  WnT2.  This  value  was  calculated  based  on  the 
assumption  that  water  flows  over  the  surface  of  the  torpedo  wall  at  a  speed  of  1m  sec'1. 

The  power  input  was  adjusted  to  get  a  maximum  temperature  of  approximately  800K 
on  the  beam.  This  was  done  in  ANSYS  by  running  the  analysis  at  a  specific  power  level 
and  obtaining  the  maximum  beam  temperatures.  Based  on  whether  the  temperature  was 
more  than  800  K  or  less  than  800  K,  the  power  to  the  beams  was  adjusted  accordingly. 
The  power  increments/decrements  were  calculated  based  on  the  relative  temperature 
difference  (+  800  -  Tmav). 

800 

Other  considerations  in  the  model: 

The  heat  is  conducted  from  the  beams  and  the  adjacent  chip  area  to  the  substrate 
through  the  2  micron  stagnant  air  gap  and  silicon-dioxide  respectively.  The  conduction 
takes  place  not  only  from  the  bottom  of  the  beam  but  also  from  the  sides  of  the  beam  to 
the  substrate.  To  account  for  the  conduction  from  the  sides,  the  conductivity  was 
multiplied  by  a  shape  factor.  The  shape  factors  for  the  air  gaps  between  the  beam  and  the 
substrate  were  calculated  from  a  separate  2-D  model.  This  included  a  rectangular  surface 
generating  heat  with  an  isothermal  surface  beneath  it  separated  by  a  2  micron  air  gap  (as 
shown  in  Figure  1 1).  The  rectangular  surface  was  surrounded  by  air  on  the  sides  also. 
The  ratio  of  the  heat  flux  from  the  sides  and  the  bottom  surface  of  the  beam  to  that  from 
just  the  bottom  of  the  beam  was  used  as  the  shape  factor. 


Fig.  11:  2  D  model  and  heat  flux  variation  to  determine  the  shape  factor  for  the 

conduction  through  the  air  gap 


This  was  repeated  for  different  air  gaps  and  the  results  are  shown  in  Figure  2-12.  The 
height  of  the  beam  was  considered  to  be  125  microns  and  the  width  was  20  microns. 
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Fig.  2-12:  Variation  of  shape  factor  with  gap  thickness  between  the  beam  and  the 

isothermal  surface 


The  model  does  not  include  radiation  effects  on  the  chip  since  these  effects  were 
found  to  be  insignificant  (of  the  order  of  0.1  K).  This  was  done  by  applying  surface  to 
surface  radiation  on  the  sides  of  the  beam  and  the  adjacent  surface  of  the  chip  (as  shown 
in  Figure  2-13),  since  this  is  where  the  maximum  temperature  difference  exists.  The 
maximum  temperature  of  the  beam  was  reduced  by  0.1  K  by  applying  radiation  in  the 
above  defined  manner. 

Natural  convection  has  not  been  included.  The  Rayleigh  number  for  this  model 
was  calculated  to  be  of  the  order  of  10'7.  The  natural  convection,  for  Rayleigh  numbers  of 
this  order,  is  very  small  and  hence  neglected. 

A  mesh  convergence  study  was  performed  for  the  S&A  chip  and  other  parts  of  the 
package  including  the  substrate  and  the  top  chip.  A  close  check  was  maintained  on  the 
beam  temperatures  while  varying  the  mesh  of  the  various  parts  of  the  package.  Based  on 
the  results,  an  optima!  mesh  was  chosen  for  the  model.  The  mesh  convergence  study 
results  can  be  seen  in  the  Figure  14. 


Surface  to  surface  radiation  applied  on  the 
sides  of  the  beam  and  adjacent  silicon 


Fig.  2-13:  Radiation  applied  to  the  surface  of  the  beam  and  the  adjacent  silicon  surface 


Mesh  convergence  -  Substrate 


Fig.  2-14:  Mesh  convergence  study  of  the  S&A  chip 


A  transient  analysis  was  performed  on  the  model.  An  automatic  time  step  was 
used  in  ANSYS  for  the  iterations.  This  means  that  the  time  step  varied  based  on  the 
performance  of  the  system.  As  the  system  started  reaching  steady  state,  the  time  step 
increased.  The  minimum  time  step  used  was  1  second.  For  the  current  simulation,  the 
areas  of  the  package  shown  below  reach  a  steady  state  in  approximately  400s  (as  shown 
in  Figure  2-15).  The  temperature  contours,  shown  later  in  the  report,  are  for  a  steady 
state  condition. 
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Fig.  2-15:  Variation  of  temperature  (different  parts  of  package)  with  time 

The  effect  of  forced  convection  is  significant  since  the  surface  area  exposed  to  the 
water  is  relatively  large.  These  effects  can  be  seen  clearly  in  Figure  2-16.  The  convection 
heat  transfer  coefficient  used  in  the  model  is  1000  Wm'2K''.  This  value  is  suitable  for  a 
water  velocity  of  approximately  1  m  sec'1.  The  temperature  of  the  water  was  assumed  to 
be  35  degrees  C  (308  K). 
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Fig.  16:  Effect  of  forced  convection  heat  transfer  coefficient  on  maximum  temperature  of 

S&A  package 


Temperature  contours  of  the  package: 

Figure  2-17  shows  the  temperature  contours  on  the  S&A  chip.  It  can  clearly  be 
seen  that  the  actuators  reach  the  desired  temperature  for  the  chosen  power.  The 
maximum  temperature  occurs  not  in  the  center  of  the  beam  but  between  the  push 
block  and  the  end  of  the  beam  as  seen  in  an  actual  powered  beam.  This  is  because  of 
the  high  heat  generation  in  the  beams  compared  to  the  push  block.  Also,  the  push 
block  has  more  surface  area  for  heat  transfer  through  the  stagnant  air  layer  relative  to 
the  heat  generation,  compared  to  the  beam.  The  maximum  temperature  as  shown  in 
the  profile  is  801.5  K. 


Fig.  2-17:  Temperature  contours  of  the  S&A  chip  including  the  actuators 

Figure  2-18  shows  the  temperature  distribution  in  the  remainder  of  the 
S&A  chip  more  clearly  since  the  actuators  were  excluded.  The  higher  temperatures, 
as  expected,  are  in  the  portion  of  the  chip  surrounding  the  arming  and  flow  sensor 
actuators,  since  the  two  actuators  are  adjacent  to  each  other  and  the  power  density  is 
higher  than  for  the  g-sensor  area. 


Fig.  2-18:  Temperature  contours  of  the  S&A  chip  excluding  the  actuators 


The  silicon-dioxide  and  the  air  beneath  the  S&A  chip  transfers  the  heat  to  the 
substrate.  The  substrate  temperatures  rise  by  approximately  50  degrees.  The  temperature 
distribution  is  similar  to  that  of  the  S&A  chip.  This  is  shown  in  Figure  2-19.  The 
maximum  temperatures  occur  right  below  the  flow  and  arming  actuators. 


The  solder  joints  near  the  flow  sensor  and  the  arming  actuators  reach  the  highest 
temperatures  (shown  in  Figure  2-20).  The  maximum  temperature  of  the  solder  joints  is 
350  K,  i.e.  42  degrees  higher  than  the  initial  temperature.  However,  this  temperature  is 
not  high  enough  to  melt  the  solder  joints  (melting  point  of  Indium,  423  K). 


Figure  2-21  shows  the  temperature  contours  of  the  top  chip  which  is  made  of 
silicon  in  the  base  case.  This  is  a  view  of  the  bottom  surface  of  the  top  chip  which  faces 
the  S&A  chip  and  is  connected  to  it  through  the  solder  joints.  The  heat  conduction  is 
therefore  mainly  through  the  solder.  The  other  form  of  heat  transfer  to  the  top  chip  is 
through  the  air  present  between  the  top  chip  and  the  device  layer  of  the  S&A  chip.  This 
distance  is  the  same  as  the  solder  height  (10  microns).  The  highest  temperatures  are 
therefore  in  the  portion  right  above  the  actuator  beams.  The  maximum  temperature  (353 
K)  is  however  less  than  that  of  the  substrate  (357  K).  This  is  because  of  the  presence  of  a 
larger  air  gap  (10  microns)  between  the  S&A  chip  and  the  top  chip  compared  to  the 
substrate  (2  microns). 


The  temperature  profile  of  the  aluminum  plate/enclosure  is  shown  in  Figure  2-23. 
The  end  where  the  S&A  package  is  located  is  the  hottest  portion.  The  forced  convection 
condition  is  applied  on  the  fuze  cap,  which  is  the  other  end  of  the  plate.  The  difference  in 
the  maximum  temperature  on  the  aluminum  plate  and  the  kovar  carrier  is  due  to  the 
presence  of  the  thermal  pads  and  board  zero,  which  offer  considerable  resistance  to  heat 
flow. 
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Fig.  2-23:  Temperature  contours  of  the  aluminum  fuze  body 


2-3.  Parametric  studies 

The  S&A  chip  requires  relatively  high  power  for  the  actuators  to  achieve  the 
necessary  displacement.  The  objective  of  the  parametric  studies  was  to  quantify  the 
power  reduction  allowed  by  various  design  changes  while  still  maintaining  a  maximum 
beam  temperature  of  800  K. 

Effect  of  silicon  dioxide  thickness: 

The  main  path  of  heat  transfer  is  from  the  beam  to  the  substrate.  The  device  layer  and 
the  handle  layer  of  the  S&A  chip  are  separated  by  a  layer  of  silicon-dioxide.  Changing 
the  thickness  of  this  layer  changes  the  thickness  of  the  stagnant  air  layer  below  the  beam. 
The  effect  of  increasing  the  thickness  of  the  oxide  layer  is  seen  in  Figure  2-24.  The  gap 
between  the  base  of  the  beam  and  the  substrate  was  varied  from  2  microns  (used  in  base 
case  model)  to  4  microns.  The  results  are  shown  as  the  power  required  per  beam  of  the 
actuator  to  reach  a  maximum  temperature  of  800K  as  a  function  of  the  gap.  An  increase 
in  the  oxide  thickness  decreased  the  power  by  a  moderate  amount.  This  decrease  is  due  to 


the  fact  that  the  main  conduction  path  for  the  heat  from  the  beam  is  down  through  the  air 
gap  to  the  substrate.  The  shape  factors  for  the  respective  gap  thicknesses  have  been 
included  in  the  calculation.  The  power  requirement  is  reduced  considerably  (.15  Watts 
per  beam)  for  an  increase  of  1  micron  in  the  gap.  However  the  effect  tends  to  reduce  as 
the  gap  is  increased  to  4  microns. 


Power  required  per  beam  of  the  actuator  vs  silicon  dioxide 

thickness 


Silicon  dioxide  thickness  (microns) 


Fig.  2-24:  Effect  of  silicon-dioxide  thickness  on  the  actuator  power  requirements 

Effect  of  top  chip  material: 

The  top  chip  in  the  base  case  was  considered  to  be  made  of  silicon.  However,  by 
changing  the  top  chip  to  glass,  the  power  required  to  reach  a  maximum  beam  temperature 
of  800  K  was  less  than  that  required  for  the  silicon  top  chip.  The  power  requirements  for 
the  three  actuators  for  both  top  chip  materials  are  shown  in  Table  2-5.  The  conductivity 
of  glass  was  considered  to  be  1  W  m'1  K1.  As  the  conductivity  of  the  top  chip  is  reduced, 
less  heat  is  transferred  from  the  S&A  chip  to  the  top  chip,  since  it  offers  a  higher 
resistance.  This  allows  the  actuators  to  reach  the  desired  temperature  for  a  smaller  power 
input.  The  glass  top  chip  reduces  the  power  requirements  per  beam  of  the  actuator  by  0.1 
Watt  compared  to  that  required  for  a  silicon  top  chip. 


Power  required  per  beam  of  the  actuator  in  Watts 

Type  of  top 
chip  material 

Flow 

Arming 

G-sensor 

Silicon 

1.95 

2.1 

1.71 

Glass 

1.81 

2 

1.65 

Table  2-5:  Power  requirements  to  reach  a  maximum  beam  temperature  of  800K  for  top 

chips  made  of  silicon  and  glass 


The  above  method  was  repeated  for  a  range  of  thermal  conductivities  to  study  the 
relative  effect  of  the  top  chip  material.  The  results  are  shown  in  Figure  2-25.  The  power 
requirements  do  not  change  considerably  for  the  thermal  conductivities  in  the  range  20  - 
150  W  m'1  K'1.  The  reduction  can  be  seen  only  in  the  range  of  1  -  20  W  m’1  K'1. 


Thermal  conductivity  of  top  chip  vs  power  requirements  of 

actuators 


Fig.  2-25:  Variation  of  actuator  power  requirements  with  top  chip  thermal  conductivity 


Effect  of  relief  cuts  in  the  top  chip: 

To  increase  the  beam’s  thermal  isolation,  relief  cuts  are  being  considered  in  the 
top  chip.  The  effect  of  these  relief  cuts  has  been  modeled  based  on  the  dimensions 
provided  by  Indian  Head  (see  Figure  2-26).  The  power  required  to  maintain  the  beam’s 
maximum  temperature  at  800K  was  determined  using  an  algorithm  to  adjust  the  power 
input  based  on  the  resulting  beam  temperature.  The  power  requirements  for  various  relief 


cut  depths  for  both  silicon  and  glass  top  chips  are  shown  in  Figure  2-27.  The  plot  shows 
the  relative  benefit  as  the  relief  cut  depth  increases. 


Fig.  2-26:  Top  chip  relief  cut  dimensions 


As  the  relief  cut  depth  in  the  top  chip  increases,  the  resistance  to  the  heat  path  through  the 
air  gap  between  the  beam  and  the  top  chip  is  increased.  The  shape  factor  increases  with 
the  air  gap,  however,  the  decrease  in  the  conductance  is  much  higher.  Hence,  less  heat  is 
transferred  from  the  beam  to  the  top  chip  which  reduces  the  power  requirements.  As  the 
cut  depth  increases  the  relative  decrease  in  the  power  reduces  for  a  given  increment. 


Effect  of  relief  cuts  in  the  substrate: 


The  relief-cut  dimensions  for  the  substrate,  provided  by  Indian  Head  (seen  in 
Figure  2-28),  were  included  in  the  model  and  the  depth  was  varied  from  10  microns  to  50 
microns.  The  base  case  model  assumptions  and  boundary  conditions  were  used  and  the 
relevant  shape  factors  to  account  for  the  conduction  from  the  sides  of  the  beam  were 
included.  The  power  input  to  the  actuator  beams  was  adjusted  such  that  a  maximum 
temperature  of  800K  was  attained  on  the  beam.  Since  the  main  path  of  heat  loss  from  the 
beam  is  through  the  air  beneath  it  to  the  substrate,  the  cut  depth  has  a  large  effect  on  the 
power  required  for  the  beam.  This  can  be  seen  in  Figure  2-29. 


Power  requirement  per  actuator  beam  to  reach  a  temperature  of 
800K  vs  relief  cut  depth  in  the  silicon  substrate 
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Fig.  2-29:  Power  requirement  per  beam  varying  with  relief  cut  depth  in  silicon  substrate 


As  the  relief  cut  depth  increases,  the  power  required  for  the  actuator  beams  to 
maintain  800K  reduces.  From  Figure  2-19  it  can  be  seen  that  after  a  relief  cut  depth  of  10 
microns,  the  reduction  in  the  power  for  a  given  depth  increment  decreases.  The  maximum 
reduction  occurs  between  no  relief  cut  and  a  relief  cut  depth  of  5  microns. 

Observations  and  Conclusions: 

The  S&A  chip  package  has  been  modeled  in  detail.  The  temperature  distribution  for 
the  whole  package  was  obtained.  Methods  to  decrease  the  power  requirements  to  the 
actuators  have  been  studied.  This  included  parametric  studies  of  potential  design  changes 
to  the  S&A  chip  package. 

Increasing  the  silicon  dioxide  thickness  resulted  in  decreasing  the  power  required  to 
maintain  the  actuator  temperature  at  800K.  Increasing  the  gap  to  3  microns  resulted  in  a 
reduction  of  0.2  Watts  per  beam  of  the  actuator  and  a  subsequent  increase  in  the  depth  by 
another  micron  resulted  in  further  reduction  of  approximately  0.1  Watt. 

The  glass  top  chip  proved  to  be  better  in  reducing  the  power  requirements  for  the 
actuators  when  compared  to  the  silicon  top  chip.  The  power  requirement  for  the  glass  top 
chip  was  0.1  Watt  less  than  the  silicon  top  chip  per  beam  of  the  actuator. 

Relief  cuts  in  the  top  chip  and  substrate  reduced  the  power  requirements  for  the 
actuators.  As  the  cut  depth  increased,  the  incremental  reduction  in  the  power  decreased. 
The  relief  cuts  in  the  substrate  had  a  larger  effect  on  the  power  reduction  compared  to 
those  on  the  top  chip  since  the  heat  conduction  path  to  the  substrate  is  a  lower  resistance 
path.  The  power  reduction  was  as  high  as  1.1  Watts  per  beam  of  the  actuator  when  the 
relief  cut  in  the  substrate  was  50  microns.  In  the  case  of  relief  cuts  in  the  top  chip,  the 
highest  reduction  in  power  was  approximately  0.3  Watt  per  beam  of  the  actuator. 

2-4  Computational  fluid  dynamics  model 
Numerical  modeling  of  the  package  -  description: 

A  computational  fluid  dynamics  model  of  the  package  was  prepared  to  determine  the  air 
pressure  and  temperature  inside  the  package.  The  CFD  model  is  similar  to  the  finite 
element  thermal  model  with  a  few  minor  changes.  The  push  blocks  were  not  included  in 
the  model  of  the  S&A  chip.  The  aluminum  block  and  torpedo  wall  were  not  included  in 
the  model,  but  are  approximated  in  a  thermally  equivalent  manner  for  a  static  analysis. 
This  was  done  by  adjusting  the  heat  transfer  coefficient  in  such  a  manner  that  it  includes 
the  resistance  to  heat  conduction  through  the  aluminum  plate  and  the  torpedo  wall.  The 
model  (shown  in  Figure  2-30)  was  built  using  FLUENT  6.1  software.  The  material 
properties  are  the  same  as  used  in  the  finite  element  thermal  model. 


CFD  Model 


Temperature  of  air  at  the  surface  of  kovar  carrier  and  the  lid: 

The  temperature  at  the  perimeter  of  the  air  volume  is  shown  in  Figure  2-31.  The 
highest  temperature  (321  K)  of  the  air  on  this  surface,  i.e.,  in  contact  with  the  lid,  occurs 
above  the  beams.  The  air  temperature  averaged  over  the  entire  cavity  in  the  package  is 
317  K. 


Fig.  2-31:  Temperature  of  the  air  at  the  perimeter  of  the  air  volume  inside  the  package. 


Temperature  of  air  inside  the  package: 


The  highest  temperature  of  air  inside  the  package  is  the  same  as  the  temperature  of 
the  beam,  since  this  air  is  adjacent  to  the  beam.  The  maximum  temperature  can  be  seen  in 
Figure  2-32.  The  average  temperature  is  calculated  by  FLUENT  using  a  volume- 
averaged  method.  The  hot  air  is  restricted  from  reaching  the  top  of  the  package,  since  the 
top  chip  absorbs  the  heat  from  the  beams.  The  conduction  path  from  the  top  chip  through 
the  solder  joints  to  the  S&A  chip  is  a  lower  resistance  pathway  when  compared  to  the 
natural  convection  on  the  surface  of  the  top  chip.  Therefore,  the  upper  surface  of  the  air 
does  not  reach  a  high  temperature  as  seen  in  the  previous  figure  (Figure  31  -  321  K). 
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Fig.  2-32:  Temperature  contours  of  air  inside  the  package  -  base  of  the  package. 


Pressure  of  air  inside  the  package: 

The  pressure  of  the  air  inside  the  package  rises  gradually  as  the  temperature  of 
air  increases  with  time.  The  pressure  of  the  air  was  approximately  3000  Pa  (0.03  atm) 
higher  than  the  initial  pressure,  when  the  system  reached  a  steady  state  condition  as  seen 
in  Figure  2-33.  The  analysis  was  performed  with  the  assumption  that  the  fluid  (air)  is  a 
compressible  ideal  gas. 
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Fig.  2-33:  Pressure  contours  of  air  inside  the  package. 


Velocity  profiles  of  air  inside  the  package: 

The  air  adjacent  to  the  beams  rises  when  the  temperature  increases  and  forms 
plumes  above  the  beams  due  to  natural  convection.  The  magnitudes  of  the  velocities  are 
small  since  the  space  above  the  beams  is  very  small.  The  maximum  velocities  are  seen 
along  the  edges  of  the  top  chip  and  the  S&A  chip.  The  air  between  the  S&A  chip  and  the 
top  chip  moves  at  a  very  low  velocity  due  to  the  restriction  in  space.  However,  as  it 
reaches  the  lateral  end  of  the  chip,  there  is  much  more  space  for  the  air  plume  to  spread. 
The  air  velocities  increase  in  this  region.  This  can  be  seen  in  Figures  34  a,  b  and  c. 

Conclusions: 

The  absolute  pressure  of  the  package  goes  up  by  approximately  3000  Pa  (0.03  atm). 
This  is  an  overall  pressure  rise  of  air  inside  the  package.  The  pressure  value  agrees  with 
that  which  would  be  calculated  by  the  ideal  gas  law  with  the  average  temperature  values 
obtained  from  the  model. 

The  average  temperature  of  the  air  inside  the  package  goes  up  approximately  9°. 
The  volume  of  air  is  large  compared  to  the  localized  heating  of  the  air  above  the  beams, 
which  results  in  a  lower  average  temperature. 


Fig.  34  a:  Velocity  profile  formed  at  g-sensor  actuator  Fig.  34  b:  Velocity  profile  formed  at  arming  actuator 


2-5  Experimental  measurements: 

The  objective  of  this  study  was  to  obtain  a  temperature  profile  of  the  chip  and  to 
validate  the  finite  element  model. 

Procedure 

The  temperature  profile  of  the  S&A  chip  was  captured  by  the  use  of  an  Infrared 
camera.  The  emissivity  of  silicon  was  given  as  a  constant  input  to  the  camera.  A 
telescopic  lens  was  used  to  capture  the  image  of  the  S&A  chip  alone. 

Calibration  of  emissivity  of  silicon 

A  silicon  chip  consisting  of  only  actuators,  excluding  the  package,  was  heated  on 
a  hot  plate  from  room  temperature  to  200°  C.  The  temperature  of  the  chip  was  monitored 


using  a  K-type  thermocouple  attached  to  the  chip.  The  emissivity  of  silicon  in  the  camera 
was  adjusted  until  the  temperature  matched  with  the  thermocouple  reading. 

Powering  the  S&A  chip 

The  wire  bonded  S&A  chip  consisting  of  two  flow  sensor  actuators,  two  arming 
actuators  and  one  g-sensor  actuator,  was  powered  using  a  DC  power  supply.  The 
package  was  not  covered  with  the  lid  to  facilitate  the  viewing  through  the  ]R  camera. 

Experimental  Setup 

The  package  was  placed  on  a  cold  plate  maintained  at  20°  C.  This  was  done  to 
maintain  an  isothermal  boundary  condition.  A  thermal  paste  was  applied  to  increase  the 
uniformity  in  the  thermal  resistance  between  the  package  and  the  cold  plate  as  seen  in  the 
Figure  2-27. 


Experimental  setup 

IR  camera 

Telescopic 
lens 

Package 

Thermal  grease 
Cold  plate 


Fig.  2-27:  Setup  for  the  experiment 


Finite  element  model 

A  few  changes  were  made  to  the  base  case  model  based  on  the  experimental 

setup. 

The  top  chip  was  not  included  in  the  model  since  the  part  did  not  have  a  top  chip. 
The  package  was  placed  on  a  cold  plate.  The  same  was  represented  in  the  thermal  model 
by  an  isothermal  boundary  condition  at  the  bottom  surface  of  the  kovar  carrier  with  a 
thermal  resistance  accounting  for  the  thermal  paste  between  the  package  and  the  cold 
plate  surface.  Natural  convection  was  included  in  the  model  on  the  walls  of  the  carrier 
and  the  surface  of  the  device  layer  of  the  chip  since  the  experiment  was  carried  out 
without  the  lid. 


The  power  in  the  chip  is  not  applied  directly  across  the  ends  of  the  actuator  beams 
in  the  experiment.  The  path  from  the  pads  to  the  beams  offers  a  resistance  to  the  current. 
This  allows  some  of  the  power  to  be  dissipated  in  this  region.  The  power  dissipated  was 
calculated  based  on  the  resistivity  values  of  the  silicon  determined  experimentally.  Since 
the  resistivity  is  temperature  dependent,  the  temperature  was  chosen  based  on  the 
approximate  temperature  reached  by  the  device  layer  of  the  chip  in  the  base  case  model. 
This  amount  of  power  was  deducted  from  the  power  being  supplied  to  the  actuator  beams 
and  was  supplied  to  the  adjacent  area. 

Averaging  technique 

The  temperature  on  the  actuator  was  captured  using  the  “point”  option  in  the  IR 
camera.  This  option  enables  the  user  to  pick  any  point  on  the  screen  and  displays  the 
temperature  of  the  chosen  location  (pixel).  The  pixel  size  for  this  case  was  50  microns. 

The  temperature  read  by  the  camera  is  an  area-averaged  temperature  over  5  pixels 
which  is  calculated  to  be  an  area  of  250  x  250  microns.  Since  the  actuator  beams  are  only 
20  microns  wide,  the  temperature  is  averaged  with  the  area  surrounding  the  silicon.  The 
actuator  beams  are  surrounded  by  air  and  hence  the  IR  camera  reads  the  temperature  of 
the  substrate  silicon  beneath  the  device  layer. 

The  same  procedure  is  adopted  to  calculate  an  area  averaged  temperature  from  the 
FEA  model.  The  temperatures  from  both  methods  for  two  different  power  levels  and  the 
respective  power  levels  are  shown  in  Table  2-6  and  Table  2-7.  The  250  x  250  micron 
area  to  be  considered  consists  of  the  actuator  beams  and  the  silicon  substrate  beneath  as 
shown  in  the  Figure  2-28. 
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Fig.  2-28:  Averaging  area  around  the  beam 
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Power  level 

Arming 

Flow 

G-sensor 

1 

1.66 

1.4 

1.08 

2 

2.325 

1.4 

.875 

Table  2-6:  Averaged  temperatures  from  the  FEA  Table  2-7:  Power  input  to  the  actuators, 

model  and  the  IR  camera 


Temperature  contours  from  the  IR  camera 

The  temperature  contours  on  the  S&A  chip  can  be  seen  in  Figure  29.  The  picture 
on  the  left  is  the  image  from  the  IR  camera.  The  temperature  shown  is  an  area  averaged 
temperature  as  explained  in  the  previous  section.  The  maximum  temperature  occurs  in 
the  portion  between  the  push  block  and  the  anchor  end  of  the  actuator  beam  similar  to  the 
contour  from  the  FEA  model.  The  push  block  is  at  a  lower  temperature  compared  to  the 
maximum  temperature.  The  picture  on  the  right  is  the  temperature  contour  from  the  FEA 
model.  However  the  temperatures  shown  are  nodal  temperatures.  As  seen  in  the  picture 
the  rest  of  the  chip  is  at  a  lower  temperature. 


Fig.  2-29:  Comparison  of  temperature  contours  from  the  IR  camera  (left)  and  the  FEA  model 

(right) 


Fig.  2-30:  Temperature  contours  excluding  the  actuators.  IR  camera  (left)  and  FEA  model 

(right) 


Figure  3-30  shows  the  temperature  contours  on  the  S&A  chip  excluding  the 
actuators.  In  the  case  of  the  1R  camera,  this  was  done  by  decreasing  the  range  of  the 
temperature  such  that  the  maximum  temperature  (occurring  on  the  beam)  falls  out  of 
range.  This  allows  for  the  thermal  gradient  to  be  seen  in  the  rest  of  the  chip.  The  solder 
joints  on  the  chip  in  the  picture  appear  to  be  colder  than  the  rest  of  the  chip.  This  is  due  to 
the  fact  that  the  emissivity  in  the  1R  camera  is  calibrated  to  that  of  silicon  and  not  indium. 
Hence  the  temperature  shown  by  the  1R  camera  for  the  solder  does  not  reflect  the  actual 
temperature. 

Conclusions 

The  following  conclusions  were  drawn  based  on  the  results  seen  from  the  IR 
measurements. 

o  Reasonably  good  agreement  has  been  obtained  between  the  experimental 
measurements  and  the  finite  element  model, 
o  The  temperatures  calculated  by  the  averaging  technique  are  quite  sensitive  to  the 
specific  area  chosen  in  the  model  and  are  therefore  subject  to  error, 
o  Good  agreement  has  been  obtained  for  the  overall  temperature  contours, 
o  These  results  give  confidence  in  the  trends  provided  by  the  parametric  studies. 

2-6  Future  work 

The  power  required  for  the  actuators  to  reach  800  K  is  still  relatively  high. 
Different  ways  of  achieving  the  purpose  of  reducing  the  power  requirement  can  be 
studied.  One  of  the  ways  to  do  this  is  to  alter  the  shape/contour  of  the  beam  to  achieve  the 
same  displacement  at  lower  temperatures  and  hence  lower  power  inputs.  Better  thermal 
isolation  of  the  beam  can  also  be  investigated. 

The  temperature  of  the  beam  is  not  known  accurately.  Hence,  the  thermo¬ 
mechanical  behavior  of  the  silicon  beams  is  difficult  to  be  verified.  Non-contact  surface 
measurements  can  be  investigated  further  to  accurately  determine  the  temperature  of  the 
actuator. 

Further  parametric  studies  can  be  carried  out  for  the  different  parts  of  the  package. 
These  parametric  studies  can  include  the  dimensions  or  the  material  properties  of  the 
package. 


Sub  Task  3:  Shock/Vibration  Durability 


Performer:  Mohammad  1  Younis,  James  Pitarresi,  Dan  Jordy 
Organization:  SUNY  Binghamton 

3-1  Dynamic  Durability  Modeling 

This  report  serves  to  document  our  progress  including  the  modeling  and 
simulation  of  the  optical  fiber,  thermal  actuator,  and  G-sensor. 

3-1-1  Optical  Fiber  -  Cantilevered  Beam 

Some  approximations  have  been  taken  during  the  analysis  of  the  optical  fibers. 
First,  ideal  fixed-free  boundary  conditions  were  assumed  to  exist,  however  in  reality  that 
may  not  be  true  [1],  Variations  in  material  properties,  due  to  fabrication  conditions,  have 
been  neglected.  The  assumed  properties  are  shown  in  Table  3-1 . 


Silica  (SiOA 

■  Length  =  4000  pm  [1] 

■  Diameter  =  125  pm  [1] 

■  Density  =  2196  kg/m3  [  1  ] 

■  Modulus  of  Elasticity  =72  GPA  [2]  [3] 

■  Poisson's  Ratio  =  0. 1 6  [2]  [3] 

■  Yield  Strength=  8.4  GPA  [2]  [3] 

Table  3-1:  Material  properties  assumed  for  the  optical  fibers.  These  properties  can  be 
expected  to  vary  because  of  process  variations. 

Lastly,  the  shock  profile  has  been  assumed  to  be  a  half-sine  shock  pulse.  This 
assumption  has  been  made  by  comparing  the  time  history  and  response  spectrum  of  both 
an  idealized  half-sine  and  saw-tooth  acceleration  shock  pulse  to  measured  data.  The 
measured  pulse  is  typical  to  those  encountered  by  dropping  an  electronic  device  on  a  hard 
surface,  as  shown  in  Fig.  3-1.  By  assuming  the  half-sine  pulse,  we  are  adding  a  factor  of 
safety  to  the  measured  data  in  the  range  of  resonance.  For  the  optical  fiber,  we  choose 
half  -  sine  pulses  of  1.0  ms  and  0.1  ms  duration  to  span  the  physical  range  for  a 
mechanical  shock. 


Fig.  3-1:  Comparison  between  idealized  shock  pulses  and  measured  data. 


Results  for  a  Clamped-Free  Microbeam 

We  determined  if  the  fibers  partially  align  under  shock  loading,  and  if  the  fibers 
break  under  shock  loads.  In  addition,  the  shock  load  necessary  to  have  the  fiber  contact 
the  substrate  was  determined.  In  order  to  validate  the  response,  several  finite  element  and 
analytical  models  were  generated.  First,  we  used  an  Euler-Bernoulli  beam  theory 
combined  with  the  shock  spectrum  of  a  spring  mass  model.  Second,  we  modeled  the  fiber 
in  ANSYS  and  applied  the  shock  as  an  applied  load.  Third,  we  used  an  equivalent  spring- 
mass  system  with  base  excitation.  Lastly,  we  used  the  same  ANSYS  model,  but  applied 
the  shock  as  a  base  excitation.  Next,  we  explain  the  models  used  and  the  results  obtained 
from  our  analyses. 


Distributed  parameter  model:  Euler  Bernoulli  beam  theory  and  shock  amplification 
factor 


A  distributed  static  load  will  displace  the  tip  of  a  cantilevered  beam  in  the 
2  a/p 

following  manner:  y  =  — - —  .  For  example,  a  600  g  distributed  static  load  will  displace 

d  E 

the  beam  5.92  pm.  For  a  dynamically  loaded 
beam,  the  maximum  displacement  will  be 
larger  than  the  static  case. 

Comparing  the  natural  period  of  the 
beam  with  the  shock  application  period,  a 
shock  amplification  factor  can  be  obtained. 

This  allows  a  static  solution  to  be  applied  to  a 
dynamic  case  [4]. 


Spectrum  shock 

t/T 


For  the  1 .0  ms  shock  pulse,  the  ratio  of  Fig.  3-2:  Shock  response  spectrum  for 
the  shock  pulse  to  the  natural  period  is  -12.5.  a  half-sine  base  acceleration.  [5] 

This  ratio  is  large  enough  so  that  we  can 
assume  a  quasistatic  response  to  the  forcing 


(amplification  factor  -1.1,  obtained  from  a  shock  response  spectrum).  The  0.1  ms  pulse 
has  an  amplification  factor  of  ~1.7.  This  factor,  found  from  Fig.  3-2,  is  multiplied  by  the 
input  force  to  account  for  the  dynamic  nature  of  the  problem. 


ANSYS  -  Distributed  load  via  structural  acceleration 

A  finite  element  simulation  program  ANSYS  was  used  to  model  the  beam. 
SOLID95  elements,  capable  of  stress  stiffening,  large  deflection,  and  large  strain  effects 
were  used.  These  elements  were  swept  through  the  length  of  the  beam  to  produce  a 
consistent  mesh.  The  acceleration  was  modeled  as  a  half-sine  pulse  with  a  duration  of 
either  1.0  ms  or  0.1  ms.  The  acceleration  was  applied  to  all  the  nodes  using  the  “ACEL” 
command  using  discrete  data  points,  solving  for  the  time  transient  solution.  Results  from 
ANSYS  are  given  in  Table  12. 


Shock  Duration  Shock  level  (g's)  Von  Mises  (MPa)  Displacement 

(micrometers) 

1-0  ms  7,452.0  126.7  125 

0.1ms  11,363.45  135.2  125 

Table  1-2:  Shock  levels  that  lead  to  a  125  pm  displacement  for  both  1.0  ms  and  0.1  ms 
shock  pulse.  Corresponding  von  Mises  stresses  are  also  shown. 

It  was  found  that  a  displacement  of  almost  8  mm  is  needed  to  fracture  the  beam, 
neglecting  contact  with  the  substrate.  This  value  is  unrealistic,  and  therefore  the  beam 
will  not  fracture  under  this  type  of  loading.  The  natural  period  for  the  beam  is  0.16  ms, 
which  is  close  to  the  shock  pulse  of  0.1  ms.  Therefore,  the  0.1  ms  shock  pulse  is  expected 
to  create  a  more  dynamic  response  in  the  fiber  than  the  1 .0  ms  shock  pulse. 

When  nonlinear  geometry  is  used  in  ANSYS,  we  observed  only  a  small  deviation 
in  displacement.  The  shock  level  investigated  for  this  deviation  was  enough  to  displace 
the  beam  -130  microns.  These  results  are  shown  in  Table  3-3. 


NONLINEAR  BEHAVIOR 


12000  g’s, 

Von  Mises  Stress 
(MPa) 

1ms 

Displacement 

(micrometers) 

Nonlinear 

130.4 

131.6 

Linear 

142.8 

132.0 

Difference  9.49% 

7500  g's  0.1ms 

0.27% 

Nonlinear 

115.4 

126.7 

Linear 

127.5 

125.8 

Difference 

10.50% 

-0.69% 

Table  3-3:  Comparison  of  the  results  for  linear  and  nonlinear  theory  for  ANSYS 
structural  acceleration.  Note  the  small  variation  in  displacement  for  a  shock  level  that 

aligns  the  fibers. 


Response  to  base  excitation 

Since  the  devices  are  very  small,  it  can  be  expected  that  they  perceive  a  shock 
pulse  as  a  base  excitation,  rather  than  as  an  applied  load.  The  last  two  models  account  for 
this  in  their  approach. 

Spring  Mass  Model 

Here  we  use  a  simple  spring  mass  model  with  an  equivalent  spring  and  mass.  The 
equation  of  motion  is  then  integrated  to  determine  its  response  to  a  base  excitation  [6]. 
Results  for  this  method  are  presented  in  0 


Comparison  between  results. 

ANSYS  modeling 

In  order  to  model  base  excitation  in  ANSYS,  a  large  mass  (m)  is  added  to  one  end 
of  the  beam.  This  mass  is  then  forced  such  that  the  structure  experiences  an  acceleration 
equivalent  to  a  =  F/m.  This  can  be  seen  in  Fig.  3-3. 


w(x,t) 


Fig.  3-3:  Schematic  representation  of 
ANSYS  base  excitation  model. 

When  nonlinear  geometry  is  used  in  ANSYS,  we  observed  only  a  small  deviation 
in  displacement.  The  shock  level  investigated  for  this  deviation  was  enough  to  displace 
the  beam  -130  microns,  which  is  the  largest  displacement  that  the  beam  can  be  expected 
to  achieve.  These  results  are  shown  in  Table  3-43-4. 


LINEAR  VS  NONLINEAR  COMPARASION 

Von  Mises  Stress  Displacement 


(MPa)  (micrometers) 

12000  g’s,  1.0  ms 

Nonlinear  (geometry)  141.2  131.7 

Linear  (geometry)  142.8  132.0 

Difference  1.17%  0.24% 

7500  g's  0.1ms 

Nonlinear  (geometry)  90.0  126.98 

Linear  (geometry)  90.2  125.98 


Difference 


.19% 


0% 


Table  3-4:  Comparison  between  linear  and  nonlinear  theory  for  ANSYS  base  excitation. 
Note  the  small  variation  in  displacement,  for  a  shock  level  that  aligns  the  fibers. 


Figure  3-5  shows  the  maximum  displacement  of  the  beam  (tip  deflection)  for  both 
the  0.1  ms  and  1.0  ms  shock  pulse.  The  0.1  ms  shock  pulse  creates  a  larger  response  in 
the  fiber  than  the  1.0  ms  pulse.  This  is  because  the  0.1  ms  shock  pulse  is  closer  to  the 
natural  period  of  the  fiber  (0.160  ms). 


Fig.  3-5:  Comparison  of  maximum  displacement  for  shock  durations  of  0.1  ms  and  1 .0 
ms.  The  red  horizontal  line  shows  the  displacement  at  which  the  optical  fibers  will  align. 


Figure  3-6  compares  the  time  history  response  of  the  beam  for  the  1.0  and  0.1  ms 
shock  pulses.  For  the  longer  shock  pulse  (1  ms),  the  beam  behaved  quasistatically  during 
the  shock  application,  and  vibrated  freely  afterward.  For  the  shorter  shock  pulse  (0.1  ms), 
the  structure  responded  as  if  it  was  subjected  to  an  impulsive  force. 


Time  history  for  1  ms  shock  pulse 


Time  history  for  0.1  ms  shock  pulse 


Fig.  3-6  :  ANSYS  results  for  the  two  shock  cases  (1.0  ms  and  0.1  ms  pulses),  applied  as  a 
base  excitation.  The  amplitude  of  the  shock  level  is  600  g’s.  Notice  the  quasi-static 
response  for  the  1  ms  shock  pulse  and  the  dynamic  response  for  the  0.1  ms  shock  pulse. 


Comparison  between  results 


Table  3-45, 

Table  3-6  and  Fig.  below  show  a  comparison  between  the  results  of  the  methods  used. 
Table  3-  shows  the  shock  level,  for  each  method  used,  to  displace  the  fiber  tip  125 
microns.  This  distance  is  enough  to  align  the  fibers  completely. 

Table  3-6  shows  the  shock  level  at  which  the  fibers  will  strike  the  substrate.  Fig. 
compares  the  shock  load  versus  tip  displacement  curves  calculated  using  all  the 
approaches. 


Amplitude  of  shock 

Ug) _ _ _ 

Model  used 

0.1  ms 
Shock 
duration 

1.0  ms 

Shock 

duration 

Euler  Beam 

7,507 

11,602 

ANSYS  Distributed  Load 

7,450 

11,363 

Spring  Mass 

8,175 

11,752 

ANSYS  Base  excitation 

7,466 

11,448 

Table  3-5:  Level  of  shock  (g’s)  to  displace  the  beam  125  microns,  aligning  the  fibers. 


Amplitude  of  shock 

(g) _ 

Model  used 

0.1  ms 
Shock 
duration 

1.0  ms 

Shock 

duration 

Euler  Beam 

KESBURI 

ANSYS  Distributed  Load 

HB  ■ 

Spring  Mass 

130.8 

ANSYS  Base  excitation 

114.4 

183.18 

Table  3-6:  Level  of  shock  (g’s)  to  displace  the  beam  2  microns,  hitting  the  substrate. 


Comparision  Between  Results  for  1.0  ms  shock  pulse 


Spring-Mass  Model  -  Base  Excitation  g  ANSYS  Base  Excitation 
Distributed  Parameter  &  Shock  Spectrum  ANSYS  -  Applied  Load 
Spring-Mass  Model,  Direct  Excitation 


a.  1.0  ms  pulse 


Comparision  Between  Results  for  1.0  ms  shock  pulse 


j  j — ♦ — Spring-Mass  Model  -  Base  Excitation  «  ANSYS  Base  Excitation  j  j 

Distributed  Parameter  &  Shock  Spectrum  ANSYS  -  Applied  Load  ! 

i  — jk — Spring-Mass  Model,  Direct  Excitation  j  j 

b.  0.1  ms  pulse 

Fig.  3-7:  A  comparison  among  various  approaches  for  the  two  shock  cases  (1 .0  ms  and 
0.1  ms  pulses).  Plotted  values  are  the  maximum  beam  displacement  for  a  given  half-sine 

shock  input. 


Parametric  Studies 


For  the  distributed  static  load  presented  in  the  first  method,  the  equation  for  the 


maximum  displacement  at  the  free  end  is 


y  = 


2a/p 

d2E 


,  where  p  is  the  density  of  the  beam 


and  ag  is  the  acceleration  applied  to  the  beam  in  units  of  length*time'2.  This  simplistic 
formula  shows  the  displacement's  dependence  on  both  the  length  and  the  diameter  of  the 
beam.  This  dependence  is  validated  with  ANSYS.  Note  that  for  the  purposes  of  this 
analysis,  linear  theory  was  used  in  ANSYS  for  the  reduced  diameter  case,  where 
significant  tip  displacement  takes  place. 


Displacement  (pm) 

7000  g’s,  1  ms 
shock  pulse 

ANSYS 

Static  Load 

Original 

76.43 

— 

1/2  length 

4.46 

4.78 

1/2  diameter 

332.30 

305.72 

3/4  diameter 

145.05 

135.87 

Table  3-7:  Design  modifications 

Deviations  from  the  expected  values  occur  because  the  natural  frequency  of  the  beam 
changes  as  the  length  and  diameter  are  varied.  This  changes  the  dynamic  nature  of  the 
ANSYS  model.  The  original  design  has  a  length  of  4000  pm  and  a  diameter  of  125  pm. 


As  the  length  is  increased  or  the  diameter  is  decreased,  the  natural  period  of  the 
structure  increases  closer  to  the  1  ms  shock  pulse.  This  creates  more  resonance  effects, 
which  makes  the  observed  displacement  in  ANSYS  higher  than  predicted  by  the  static 
case,  where  this  phenomenon  is  not  accounted  for.  The  opposite  effect  can  be  seen  when 
halving  the  length.  As  the  natural  frequency  gets  closer  to  the  shock  pulse,  the  shock 
amplification  factor  increases,  which  further  creates  variations  between  the  predicted 

static  case  and  the  dynamic  case.  Therefore,  the  general  trend  of  y  00  and  yK/4  from 
the  static  case  can  only  be  used  as  a  first  approximation. 


Effect  of  damping 

Damping  in  the  system  can  be  modeled,  in  the  spring  mass  system,  by  including  a 
term  that  is  proportional  to  the  velocity  of  the  mass.  In  Table  and  Table  we  show  the  tip 
displacement  that  occurs  for  varying  levels  of  damping.  For  the  highest  damping  ratio,  it 
can  be  seen  that  the  tip  displacement  for  the  0.1  ms  shock  pulse  is  lower  than  that  of  the 
1 .0  ms  pulse.  This  is  opposite  from  what  is  seen  when  damping  is  not  present.  As  the 
damping  ratio  increases,  the  shock  response  spectrum  decreases.  At  a  damping  ratio  of 
0.5,  the  amplification  factor  is  lower  in  the  0.1  ms  case  than  in  the  1.0  ms  case. 


Maximum  Tip  Deflection 
(micrometers)  for  a  shock  level  of 
1500  g’s 

1.0  ms 

0.1  ms 

Damping  Ratio 

0 

16.11 

23.15 

15.32 

21.45 

0.5 

14.76 

12.76 

Table  3-8:  Effect  of  damping  in  a  spring  mass  model  with  base  excitation  for  1 .0  and  0.1 
ms  shock  pulses,  at  a  shock  level  of  1 500  g’s. 

Notice  that  for  heavy  damping,  the  response  to  the  0.1  ms  shock  duration  is  lower  than 
the  response  to  the  1 .0  ms  shock  pulse. 


Maximum  Tip  Deflection 
(micrometers)  for  a  shock  level  of 
5000  g’s 

1.0  ms 

0.1  ms 

Damping  Ratio 

0 

53.71 

77.18 

■BB 

51.06 

71.49 

0.5 

49.21 

42.53 

Table  3-9:  Effect  of  damping  in  a  spring  mass  model  with  base  excitation  for  1 .0  and  0.1 
ms  shock  pulses,  at  a  shock  level  of  5000  g’s. 

Notice  that  for  heavy  damping,  the  response  to  the  0. 1  ms  shock  duration  is  lower  than 
the  response  to  the  1 .0  ms  shock  pulse. 


Mechanical  Stress 

Optical  fibers  produced  from  synthetic  fused  silica  have  a  theoretical  strength  of 
~2,000  kpsi  (-13.8  GPa),  based  upon  the  Si-O  bond  strength.  The  presence  of 
imperfections  and  flaws  in  the  bulk  of  the  material  and  on  its  surface  lowers  the  observed 
strength  to  ~700  kpsi  (4.8  GPa).  Other  literature  states  the  fracture  strength  of  silica  to  be 
8.4  GPa  [3]. 

Table  1  shows  that  the  von  Mises  equivalent  stress,  at  a  level  capable  of  aligning 
the  fiber,  would  not  cause  the  fiber  to  fracture. 


3-1-2  Thermal  Actuator  Modeling 

Some  approximations  have  been  taken  during  the  analysis  of  the  thermal  actuator. 
First,  ideal  clamped-clamped  boundary  conditions  were  assumed,  however  in  reality  this 
may  not  be  true  [1],  Variations  in  material  properties,  due  to  fabrication  conditions,  have 
been  neglected.  In  addition,  single  crystal  silicon  behaves  anisotropically,  and  it  is 
therefore  difficult  to  obtain  single  values  for  the  material  properties  that  are  used  with 


certain  models.  Literature  gives  a  range  of  Young's  Moduli,  including  190  GPa  [9],  166 
GPa  [1],  and  150  GPa  (u=0.17)  [10].  The  properties  assumed  throughout  this  analysis  are 
shown  in  Table . 


Silicon  (Si) 

■  Length  =  5500  pm  [1]  [8] 

■  Thickness  =  125  pm  [1]  [8] 

■  Width  =  20  pm  [1]  [8] 

■  Offset  =  45  pm  [1]  [8] 

■  Density  =  2332  kg/mJ[l]  [9] 

■  Modulus  of  Elasticity  =  1 54  GPa 

■  Poisson's  Ratio  =  0.3 

■  Yield  Strength=  7  GPA  [9]  [3] 


Table  3-10:  Material  properties  assumed  for  the  thermal  actuator.  These  properties  can 
be  expected  to  vary  because  of  process  and  material  variations. 

Again,  a  half-sine  shock  pulse  has  been  used  to  represent  the  shock  profile  (See 
Fig.  3-1  and  Section  Error!  Reference  source  not  found.).  We  choose  sine  pulses  of  1.0 
ms  and  0.1  ms  duration  to  span  the  physical  range  for  a  mechanical  shock. 


Results  for  a  Clamped-Clamped  Thermal  Actuator 

In  order  to  determine  its  response  to  shock  pulses  of  1.0  and  0.1  ms  and  to 
minimize  computational  requirements,  we  modeled  the  thermal  actuator  with  half¬ 
symmetry.  We  verify  if  the  thermal  actuator  is  capable  of  pushing  against  the  optical 
fiber  during  shock  loads.  In  addition,  we  establish  the  shock  load  necessary  to  have  the 
thermal  actuator  contact  the  substrate. 

Finite-Element  Model 

We  used  ANSYS,  a  finite  element  simulation  program,  to  simulate  the  device’s 
response  to  shock  loads.  SOLID95  elements  were  used  to  model  the  thermal  actuator. 
These  elements  are  capable  of  capturing  nonlinear  behaviors,  such  as  stress  stiffening, 
large  deflection,  and  large  strain  effects.  The  elements  were  swept  through  segments  of 
the  thermal  actuator  to  produce  a  consistent  mesh.  The  thermal  actuator  was  modeled  as  a 
clamped-clamped  frame  with  an  offset,  with  a  center  section  containing  regular 
perforations  (25pm  x  25pm,  with  25  pm  spacing  between  each  perforation).  The  shock 
was  modeled  as  a  half-sine  pulse  with  a  duration  of  either  1.0  ms  or  0.1  ms.  The  shock 
was  applied  to  all  the  nodes  of  the  finite  element  model  via  the  “ACEL”  command  using 
discrete  data  points,  solving  for  the  time  transient  solution.  In  our  first  report  on  the 
modeling  and  simulation  of  the  optical  fiber,  we  show  that  applying  a  distributed  load 
across  a  structure  is  comparable  to  a  base  excitation  on  the  structure.  Since  the  small 
structures  are  expected  to  perceive  a  shock  pulse  as  a  base  excitation,  and  we  have 
previously  shown  that  a  distributed  load  across  a  structure  is  equivalent  to  a  base 


excitation,  then  we  are  justified  in  using  such  approach  here.  Fig.  through  Fig.  show  the 
beam  without  symmetry,  and  with  the  length,  offset,  thickness  and  width  labeled. 


Fig.  3-9:  Closer  view  of  offset 


In-Plane  Results 


We  study  the  response  of  the  thermal  actuator  to  in-plane  shock  loads  of  duration 
1.0  ms  and  0.1  ms.  Fig.  shows  a  comparison  between  the  results  of  linear  and  nonlinear 
theories,  revealing  that  geometric  nonlinearities  are  significant  even  for  small  shock 
levels.  We  conclude  that  a  nonlinear  beam  theory  and  analysis  has  to  be  used  to  simulate 
the  mechanical  response  of  the  actuator  to  in-plane  loads.  It  is  clear  from  Fig.  that  the 
response  due  to  the  0.1  ms  shock  load  is  larger  than  that  of  the  1.0  ms  shock. 
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Fig.  3-11:  Comparison  between  linear  and  nonlinear  theory  for  ANSYS  distributed  load. 
Note  the  variation  in  displacement  between  linear  theory  and  nonlinear  theory  at 

increasing  shock  levels. 


To  justify  this,  we  calculate  the  natural  frequency  of  the  actuator  (Table).  The  natural 
period  for  the  thermal  actuator  in  the  in-plane  direction  is  0.094  ms,  which  is  close  to  the 
shock  pulse  of  0.1  ms.  Therefore,  the  actuator  exhibits  a  dynamic  response  to  a  0.1  ms 
shock  pulse,  and  displaces  further  than  when  excited  by  the  longer  shock  pulse  of  1 .0  ms. 


Description  of  mode 

0.0939 

In-plane  motion  -  Looks  similar  to  the  response  to  an 
in-plane  load 

0.0605 

Twisting  about  the  actuator  arms  (clamped-clamped 
beams)  -  Looks  similar  to  the  response  to  an  out-of¬ 
plane  load. 

0.0445 

Actuator  arm  stationary,  stacked  clamped-clamped 
beams  moving  in  opposite  directions  in-plane 

0.0362 

In-plane  motion  -  Looks  similar  to  the  3rd  mode  of  a 
clamped-clamped  beam 

0.0324 

Out-of-plane  motion  -  Twists  about  middle  of 
actuator  arm  -  Center  of  clamped-clamped  beams 
displaces  out-of-plane 

Table  3-11:  First  5  natural  periods  for  the  thermal  actuator. 


a)  First  Mode  Shape 


b)  Second  Mode  Shape 


c)  Third  Mode  Shape  d)  Fourth  Mode  Shape 


Fig.  3-12:  First  5  mode  shapes  of  the  thermal  actuator 


In  Fig.  ,  we  compare  the  response  of  the  thermal  actuator  to  that  of  the 
cantilevered  optical  fiber  under  the  same  shock  load.  We  notice  that  a  significant  shock 
load  is  necessary  to  force  the  actuator  in  the  in-plane  direction,  compared  to  the  optical 
fiber.  Therefore,  we  conclude  that  a  shock  load  in  the  in-plane  direction,  towards  fiber 
alignment,  will  not  cause  the  thermal  actuator  to  contact  the  optical  fiber. 


Fig.  3-13:  Response  of  optical  fiber  and  thermal  actuator  to  in-plane  shock  loads  of  0. 1  & 
1 .0  ms  durations.  The  thermal  actuators  will  be  incapable  of  pushing  the  optical  fibers 

into  alignment  during  shock  loads. 


Out-of-Plane  Results 

We  investigate  the  response  of  the  thermal  actuator  to  an  out-of-plane  shock  load. 
Figure  3-4  shows  the  response  of  the  thermal  actuator  to  shock  loads  of  durations  0.1  ms 
and  1 .0  ms,  calculated  using  linear  and  nonlinear  theories.  We  observe  that  in  the  out-of¬ 
plane  direction,  a  displacement  of  4  pm  does  not  cause  significant  nonlinearities  to  arise. 
Hence  the  thermal  actuator  behaves  linearly  in  this  direction.  For  a  1.0  ms  shock 
duration,  a  shock  level  of  approximately  1,500  g’s  is  needed  to  cause  the  device  to  hit  the 
substrate,  based  on  a  2  micrometer  gap,  and  approximately  3000  g’s  to  hit  the  substrate 
based  on  a  4  micrometer  gap.  For  the  0.1  ms  shock  duration,  a  shock  level  of 
approximately  1200  g’s  is  needed  to  cause  the  device  to  hit  the  substrate,  based  on  2 
micrometer  gap,  and  approximately  2250  g’s  to  hit  the  substrate  based  on  a  4  micrometer 
gap.  Figure  3-4  shows  that  the  response  to  the  shorter  shock  pulse  (0.1  ms)  is  larger  than 
the  response  to  the  longer  shock  pulse  (1.0  ms).  This  can  be  justified  by  observing  that 
the  second  natural  period  (Table)  is  close  to  the  0.1  ms  shock  duration,  and  that  the 
second  mode  shape  is  similar  to  the  response  to  an  out  of  plane  shock.  Similar  to  the  in¬ 
plane  case,  the  actuator  exhibits  a  dynamic  response  to  the  0.1  ms  shock  pulse,  and 
displaces  further  than  when  excited  by  the  1 .0  ms  shock  pulse. 

In  addition,  we  find  that  the  tip  of  the  thermal  actuator  will  contact  the  substrate 
first.  Fig.  shows  the  deformed  shape  of  the  actuator  when  subjected  to  an  out  of  plane 
shock  of  2000  g’s,  for  a  shock  duration  of  0.1  ms.  The  maximum  displacement  in  the  out 
of  plane  direction  is  approximately  3.6  micrometers  in  this  case.  This  can  be  justified  by 
visualizing  the  thermal  actuator  as  a  clamped-clamped  beam,  with  a  cantilevered  beam 
attached  to  its  center  and  extending  perpendicularly  from  it.  When  both  structures  are 


subjected  to  a  distributed  load,  the  clamped-clamped  beam  will  deflect  the  greatest 
distance  at  its  center,  where  the  cantilevered  beam  is  located.  The  tip  of  the  cantilevered 
beam  will  deflect  even  more.  This  visualization  justifies  the  observation  that  the  tip  of  the 
thermal  actuator  touches  the  substrate  first. 


Figure  3-4:  Out-of-plane  displacement  due  to  shock  load. 


AIMSYS 


Fig.  3-15  Deformed  shape  due  to  out  of  plane  shock  of  2000g.  Also  shown  is  the  outline 
of  the  unforced  actuator. 

Parametric  Studies 

Several  design  parameters  have  been  identified  that  are  capable  of  increasing  the 
thermal  actuator’s  resistance  to  shock.  For  an  in-plane  shock  load  these  include 


increasing  the  width  of  the  thermal  actuator  beam,  and  decreasing  the  length  of  the  beam. 
For  an  out-of-plane  shock  load  these  include  increasing  the  thickness  of  the  thermal 
actuator  beam  and  decreasing  the  length  of  the  beams,  as  defined  in  Fig.  .  Also 
investigated  was  the  number  of  clamped-clamped  thermal  arms  attached  to  the  center 
actuator  arm.  The  offset  in  the  center  of  the  beams  has  been  found  to  influence  the 
resistance  to  shock. 

Fig.  shows  that  decreasing  the  length  increases  the  in-plane  shock  resistance  of 
the  microbeam.  Fig.  shows  the  response  of  the  beam  to  out-of-plane  shock  loads  when 
the  length  of  the  thermal  actuator  is  varied. 
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Fig.  3-16:  Effect  of  varying  the  length  of  the  thermal  actuator  on  the  shock  resistance  of 
the  thermal  actuator  to  in-plane  shock  loads. 


Fig.  3-17:  Effect  of  varying  the  length  of  the  thermal  actuator  on  the  shock  resistance  of 
the  thermal  actuator  to  out-of-plane  shock  loads. 


The  addition  of  actuator  beams  (stacking)  will  decrease  the  response  of  the 
structure  to  a  given  shock  load.  At  a  shock  level  of  7000  g's  and  0.1  ms  shock  duration,  a 
thermal  actuator  with  a  single  clamped-clamped  beam  will  have  a  displacement  of  24 
pm.  For  the  same  shock  load,  but  with  two  clamped-clamped  beams,  the  maximum 
displacement  will  be  19  pm.  Therefore,  by  increasing  the  number  of  clamped-clamped 
beams,  we  can  increase  the  robustness  of  the  device  to  shock  and  vibration.  Fig.  shows 
this  effect.  It  also  shows  the  effect  of  varying  the  length  of  the  beam.  Additional 
robustness  could  be  obtained  by  implementing  a  “cascading"  thermal  actuator  device 
[13],  however  this  would  necessitate  a  major  redesign  of  the  chip,  and  no  work  has  been 
done  analyzing  this  type  of  arraignment. 


80 

70 

60 

50 

40 

30 

20  ■ 

10 


- 1  Arm 

2  Arms  (As  designed) 
—  -  -3  Arms 
4  Arms 


/ 

/ 

/ 

/ 

/ 


/ 

/ 

/ 

/  ■/ 


2000  3000  4000  5000  6000  7000  8000  9000 

Thermal  Actuator  Length  (micrometers) 


Fig.  3-18:  Effect  of  the  addition  of  actuator  beams  and  their  length  on  the  maximum 
response  to  an  in-plane  shock  load  of  1000  g's  and  0.1  ms  duration 


Fig.  shows  the  effect  of  the  offset  on  the  response  of  the  beam.  Here  the  duration 
of  the  shock  load  is  0.1  ms.  This  figure  illustrates  that  a  larger  offset  increases  the 
device’s  resistance  to  shock  in  the  in-plane  and  out-of-plane  directions.  Previous  papers 
[12]  have  shown  that  increasing  the  offset  will  in  turn  decrease  the  tip  displacement 
during  thermal  actuation.  These  two  effects  need  to  be  balanced,  such  that  proper  thermal 
actuation  can  take  place  while  ensuring  vibration  and  shock  robustness. 


beam  to  an  out  of  plane  distributed  static  load,  where  the  displacement  is  ymm  =  » 

where  1  =  ~(width)(thickness)\  A  shock  amplitude  of  15,000  g’s  is  needed  to  displace 

a  145  pm  thick  actuator  the  same  distance  as  a  125  pm  at  10,000  g’s.  Therefore, 
increasing  the  thickness  25  pm  from  the  original  design  will  yield  an  improvement  of 
50%  in  shock  resistance.  Note  that  this  change  may  not  be  feasible,  as  it  would  entail 
etching  the  features  of  the  device  further  into  the  substrate,  which  would  change  other 
features  on  the  chip.  We  have  not  focused  on  increasing  the  shock  resistance  in  the  in¬ 
plane  direction  by  increasing  the  width,  as  the  thermal  actuator  is  sufficiently  stiff  enough 
not  to  contact  the  optical  fiber  during  shock  loading. 


Thermal  Actuator  Thickness  {out  of  plane  direction)  (micrometers) 


Fig.  3-20:  Effect  of  varying  the  thickness  in  the  thermal  actuator. 

Mechanical  Stress 

The  fracture  strength  of  silicon  is  7  GPa  [9]  [3].  In  order  to  examine  the  stress 
profile  in  the  thermal  actuator,  we  displaced  a  node  on  the  tip  of  the  beam  125  pm  with 
ANSYS.  In  doing  so,  we  find  a  maximum  von  Mises  stress  in  the  actuator  of  2.3  GPa. 
Furthermore,  this  stress  occurred  at  the  point  of  application  of  the  displacement,  and  the 
stress  around  this  point  dropped  away  quickly.  Since  the  observed  stress  is  much  lower 
than  the  fracture  strength  of  the  silicon,  and  may  have  been  artificially  increased  by  the 
presence  of  the  point  displacement,  we  determine  that  the  thermal  actuator  will  not 
fracture  at  this  displacement. 


3-1-3.  G-sensor 


Here  we  determine  the  effect  that  the  squeezed-film  damping  effect  on  G-sensor 
during  out  of  plane  motion,  and  determine  the  effect  of  design  parameters  on  the  squeeze 
film  damping  effect. 


Structural  Modeling 

In  order  to  model  the  G-sensor’s  response  to  shock  and  squeeze  film  damping,  we 
first  separated  the  structural  and  fluidic  domains  to  analyze  each  separately.  For  the 
structural  domain,  we  modeled  the  G-Sensor  in  ANSYS  using  4-Node  Finite  Strain  Shell 
elements  (SHELL181).  This  element  is  suitable  for  moderately  thick  shell  structures,  and 
is  able  to  model  linear,  large  rotation,  and  large  strain  nonlinear  problems.  This  element 
is  also  capable  of  displacements  in  3  directions  at  each  node,  making  it  suitable  for  out  of 
plane  displacements.  Fig.  shows  the  model,  which  includes  the  holes  and  the  springs. 
The  G-sensor  is  constructed  out  of  silicon,  whose  properties  are  listed  in  Table  .  This 
model  does  not  have  any  restrictions  on  its  movement;  the  G-sensor  is  a  flexible  plate 
that  can  deform  in  any  direction. 


Fig.  3-21:  Model  of  the  G-Sensor,  structural  domain. 

The  red  dots  indicate  points  where  the  displacement  is  measured  in  Fig.  and  Fig. . 

Modal  Analysis 

We  first  performed  a  modal  analysis  in  ANSYS  on  the  G-sensor,  to  determine  its 
mode  shapes  and  natural  frequencies.  Fig.  shows  the  first  4  out-of-plane  mode,  and 
Table  shows  their  natural  periods.  This  analysis  illustrates  that  there  are  several  natural 
frequencies  within  a  close  range,  so  it  is  possible  that  several  modes  will  be  excited  as  the 
G-sensor  responds  to  shock  loading. 


Is' mode:  2075.9  Hz -0.48  ms  2nd  Mode:  2381.1  Hz  -0.42  ms 


3rd  Mode:  3248.4  Hz-0.31  ms  4,h  Mode:  24807  Hz- 0.04  ms 

Fig.  3-24:  First  4  out-of-plane  mode  shapes  and  their  corresponding  natural  frequencies 


Mode 

Natural 

Frequency 

(Hz) 

Natural 

Period 

(ms) 

Shock  ratio  Tpuise/Tna, 

1.0  ms 
shock 

0.1  ms 
shock 

0.5  ms 
shock 

1 

2075.9 

0.48 

4.17 

0.42 

2.08 

2 

2381.1 

0.42 

4.76 

0.48 

2.38 

3 

3248.4 

0.31 

6.45 

0.65 

3.23 

4 

24807 

0.04 

50 

5 

25 

Table  3-12:  First  4  out-of-plane  natural  frequencies  and  periods. 


Mode 

Natural 

Frequency 

(Hz) 

Natural 

Period 

(ms) 

Shock  ratio  Tpuise/Tnal 

1.0  ms 
shock 

0.1  ms 
shock 

0.5  ms 
shock 

1 

214.16 

4.669 

0.42832 

0.042832 

0.21416 

2 

1993.3 

3.9866 

0.39866 

1.9933 

3 

2133.4 

4.2668 

0.42668 

2.1334 

4 

esezih 

4.9546 

0.49546 

5 

refill 

0.404 

4.9562 

6 

4.9572 

7 

2489 

0.402 

4.978 

8 

6970.8 

0.143 

13.9416 

1.39416 

6.9708 

Table  3-13:  First  8  non-  out-of-plane  natural  frequency  and  periods. 


Undamped  Shock  Response 

In  order  to  characterize  the  response  of  the  G-sensor  to  shock  loading,  we  first 
analyze  the  device  without  damping.  By  neglecting  all  damping,  we  are  assuming  a 
worst-case  scenario.  This  also  provides  a  baseline  to  compare  to  the  damped  model.  As  in 
the  previous  work,  we  modeled  the  shock  as  a  half-sine  pulse,  applying  it  at  discrete  time 
intervals  by  using  the  “ACEL”  command  in  ANSYS.  See  Fig.  3-1  and  Section  0  for  an 
explanation  of  this  assumption. 


10g,  half  sine  shock  duration 


Fig.  3-23:  Maximum  displacement  of  the  G-sensor  to  a  shock  pulse  of  lOg  and  varying 
duration  (left).  The  maximum  point  corresponds  to  the  peak  of  the  shock  spectrum 

(right). 


As  before,  we  needed  to  choose  a  shock  duration  that  created  a  large  displacement 
in  the  G-Sensor.  We  found  this  duration  to  be  approximately  0.3  ms.  Fig.  shows  these 
results.  Fig.  shows  the  shock  amplitude  needed  to  cause  the  G-sensor  to  contact  the 
substrate  for  different  gap  widths.  It  is  clear  that  without  damping,  a  very  low  shock  level 
is  capable  of  causing  the  G-Sensor  to  contact  the  substrate. 


Shock  threshold  to  hit  the  substrate 
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Fig.  3-24:  Shock  amplitude  needed  to  cause  contact  with  the  substrate,  for  both  a  2  pm 

and  4  pm  gap. 

Fig.  and  Fig.  show  the  response  at  each  comer  of  the  G-Sensor  to  an  out-of- 
plane  shock  of  10  g’s.  Here  the  period  of  the  smaller  oscillations  is  about  0.5ms,  which  is 
close  to  period  of  the  first  3  out-of-plane  modes  (Table  ).  The  larger  period  seen  in  Fig. 
is  due  to  the  constructive  interference  between  the  first  two  out  of  plane  modes.  The 
difference  between  the  first  two  natural  frequencies  is  approximately  200  Hz,  which  is 
the  same  as  the  outer  envelope  of  the  response. 


Response  of  G-Sensor  to  a  10-g  0.1  ms  out-of-plane  shock 


Fig.  3-25:  Response  of  the  G-Sensor  to  an  out  of  plane  shock  of  10  g's  at  a  duration  of 
0.1  ms.  The  maximum  displacement  is  0.86  pm.  See  Fig.  for  the  physical  locations  of  the 

plotted  points  on  the  G-Sensor. 


Response  of  G-Sensor  to  a  10-g  1  0  ms  out-of-plane  shock 


Fig.  3-26:  Response  of  the  G-Sensor  to  an  out  of  plane  shock  of  1 0  g’s  at  a  duration  of 
1 .0  ms.  The  maximum  displacement  is  0.85  pm.  See  Fig.  for  the  physical  locations  of  the 

plotted  points  on  the  G-Sensor. 


Modeling  of  Squeeze  Film  Damping 
Theory 


At  the  micrometer  length  scale,  the 
movement  of  a  proof  mass  normal  to  a 
stationary  wall  squeezes  the  air  gap  between 
the  two  objects.  This  causes  significant  forces 
from  the  air  gap  on  the  proof  mass,  affecting 
its  dynamics.  This  force  can  have  a 
component  in  phase  with  the  displacement, 
and  a  component  in  phase  with  the  velocity. 
The  former  force  leads  to  a  restoring  spring 
force,  while  the  latter  creates  a  damping  force. 
At  low  frequency  movements,  the  air  is  able  to 
move  out  from  under  the  structure, 
minimizing  the  spring  force  and  allowing  the 

damping  to  become  the  dominant  force  from 
the  air  gap.  At  higher  frequencies,  the  air  is 
unable  to  move  out  from  under  the  device, 
creating  a  spring  force  and  lowering  the 
damping.  [14]  This  effect  can  be  seen  in  Fig.  . 


Fig.  3-27:  Damping  and  stiffness  coefficients 
from  squeeze  film  damping  for  varying 
frequencies. 


The  shape  of  the  structure  influences  the  effect  of  squeeze  film  damping  by 
affecting  the  venting  of  the  gas  from  underneath  the  device.  The  addition  of  perforations 
into  the  device  also  affects  the  squeeze  film  damping  effect  by  allowing  air  to  move  out 
through  the  perforations. 

Since  the  G-Sensor  is  thick  and  supported  by  compliant  springs,  we  can  model  the 
damping  effect  by  assuming  the  accelerometer  to  be  a  rigid  plate  moving  normal  to  the 
substrate.  Here  the  G-Sensor  can  only  move 
up  and  down,  with  no  tipping  motion 
allowed.  This  assumption  allows  us  to 
decouple  the  structural  and  fluidic  domains. 

We  modeled  the  thin  fluid  domain  using  3-D 
squeeze  film  fluid  elements  (FLU1D136). 

This  element  models  the  flow  of  a  viscous 
fluid  in  small  gaps  that  are  bound  by  a  fixed 
surface  and  a  surface  moving  perpendicularly 
to  the  fixed  surface.  This  element  is  based  on 
the  Reynolds'  Equation,  and  assumes  isothermal  behavior,  small  pressure  changes 
compared  to  the  ambient  pressure,  and  small  displacements  compared  to  the  film 
thickness.  3-D  viscous  fluid  link  elements  (FLUID]  38)  were  used  to  model  the  viscous 
flow  through  the  perforations  in  the  device. 


FLUID138 


Fig.  3-28:  Schematic  used  to  model 
squeeze  film  damping. 


Boundary  Conditions  -  Perforations 


No  Holes 


There  are  three  different  ways  in  which  _ 

ANSYS  can  model  the  perforations  and  their  specified  y 

effect  on  squeeze  film  damping.  First,  the  /§' 

perforations  can  be  neglected  if  they  are  FLUID136 

either  very  long  or  very  slender.  Physically,  Fig.  3-29:  No  holes/infinite  resistance 
this  means  that  there  is  a  very  large  flow  pressure  boundary  condition  for  the  holes, 
resistance  through  the  hole,  and  consequently 

the  pressure  underneath  each  perforation  is  approximately  the  same  value  as  for  the  un¬ 
perforated  case.  This  case  is  equivalent  to  having  no  holes  in  the  G-Sensor.  Fig.  shows 
this  case. 


No  Hole  Resistance 

If  the  perforations  are  very  thin  or  , 

.......  .  ,  Zero  Pressure  BC 

very  wide,  there  is  little  resistance  to  the  ^ ^ 

flow  through  a  hole.  Here,  the  fluid 

underneath  the  perforation  has  a  pressure  5 

equal  to  the  ambient  pressure.  To  Fig.  3.30:  No  hole  resistance  pres 

accomplish  this  in  ANSYS,  we  set  the  boundary  condition  for  the  holes. 

pressure  boundary  condition  along  every 

hole  to  be  equal  to  the  ambient  pressure.  Fig.  illustrates  this  case. 


FLUID136 


Fig.  3-30:  No  hole  resistance  pressure 


Finite  Fluid  Resistance 


A  third  way  to  model  the  squeezed  Zero  PressureBC 
film  effect  is  to  include  the  finite  H§|j  FLUID138 

resistance  in  the  perforations.  This  is  the 

most  accurate  way  to  model  the  FLU1D136 

perforations,  as  it  accounts  for  the  pressure 

drop  through  each  perforation.  To  do  this.  Fig.  3-31:  Finite  Pressure  boundary 
we  introduce  a  new  element,  FLUID138,  condition  for  the  holes, 
which  models  the  flow  through  one 

perforation.  Here,  the  pressure  at  one  end  of  the  FLUID!  38  element  is  coupled  to  the 
surrounding  FLUID136  elements,  and  at  the  other  end  the  pressure  is  set  to  the  ambient 
pressure.  Fig.  demonstrates  this  case. 


Compressible  Fluid  Assumption 

In  order  to  account  for  the  compressibility  of  the  gas,  we  subject  the  fluid  to  a 
harmonic  fluid  velocity.  By  varying  the  frequency,  we  can  examine  how  the  fluid  will 
react  when  the  G-sensor  vibrates.  We  notice  that  the  squeeze  film  damping  coefficient 
does  not  vary  greatly  from  0  Hz  to  50  000  Hz,  while  the  spring  force  increases  slightly 
from  zero.  Since  the  damping  force  is  greater  than  the  spring  force  over  this  frequency 
range,  the  fluidic  damping  can  be  considered  constant,  and  the  compressibility  of  the 
fluid  can  be  neglected  [15], 


Incompressible  fluid  assumption 

By  neglecting  the  compressibility  of  the  fluid,  we  simplify  the  problem  by 
neglecting  the  spring  effect  from  the  fluid  and  considering  the  damping  effect  only.  To  do 
end,  we  subject  the  fluid  elements  to  an  arbitrary  velocity,  which  creates  a  pressure 
distribution  in  the  fluid.  Integrating  the  pressure  distribution  over  the  area  determines  the 

F 

damping  force,  which  can  be  used  to  determine  a  damping  coefficient  by  c  =  —  [15]. 

v. 

This  damping  coefficient  is  the  c  value  in  the  single  degree-of-freedom  equation 

c 

mx  +  cx  +  kx  =  /(/) .  It  can  be  converted  to  a  damping  ratio  by  the  formula,  C,  = - , 

2  mcon 

where  con  is  the  first  out-of-plane  natural  frequency  in  rad/s.  This  value  is  a  first 
approximation,  since  the  real  structure  will  not  oscillate  at  its  first  out  of  plane  natural 
frequency,  like  a  one  degree-of-freedom  system.  In  this  application,  a  low  damping  ratio 
is  not  desirable,  as  it  may  allow  the  g-sensor  to  contact  the  substrate  during  shock 
loading,  which  could  damage  the  sensor  or  render  it  inoperable  due  to  stiction  effects. 


Hole  Boundary  Condition 

Damping  Ratio 

No  Holes 

309 

No  Hole  Resistance 

0.15 

Finite  Resistance,  25  pm  perforations 

0.21 

Table  3-14:  Damping  Ratios  for  different  boundary  conditions. 


Results 


For  the  first  case,  where  the  holes  are  neglected,  the  damping  ratio  is  found  to  be 
309.  Fig.  shows  the  pressure  distribution  for  this  case,  which  is  highest  in  the  center  and 
decreases  to  ambient  on  the  periphery. 

In  the  case  of  no  hole  resistance,  the  damping  ratio  decreases  to  0.1453.  Fig. 
shows  the  pressure  distribution  in  the  G-Sensor  for  the  case  of  no  hole  resistance.  The 
holes  can  be  seen  to  lower  the  pressure  exerted  on  the  G-sensor,  thereby  decreasing  the 
damping  ratio.  Including  the  hole  resistance  increases  the  damping  ratio  slightly  to 
0.2125. 
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Fig.  3-32:  Pressure  distribution  due  to  squeeze  film  damping  -  no  boundary  conditions  on 
the  perforations.  A  gap  thickness  of  2  pm  is  assumed. 
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Fig.  3-33:  Pressure  distribution  due  to  squeeze  film  damping  -  ambient  pressure  on  the 
perforations.  A  gap  thickness  of  2  pm  is  assumed.  Notice  the  lower  pressures  compared 

to  Fig. . 


Effect  of  hole  size 


The  high  density  of  perforations  allows 
the  fluid  to  escape  from  underneath  the  G- 
sensor  easily,  creating  a  low  damping  force.  In 
order  to  increase  the  damping  in  the  structure  to 
prevent  stiction,  we  suggest  that  the  area  of  the 
perforations  in  the  device  be  minimized.  A  hole 
size  of  10  micrometers  leads  to  a  damping  ratio 
of  4.2712,  which  is  a  significant  increase  from 
0.2125.  There  are  limits  to  this,  as  the 
perforations  appear  to  be  needed  in  order  to 
release  the  structure  from  the  substrate  during 
fabrication. 


_ § _ 

Gap 

Width 

2  pm 

4  pm 

25  um 

0.2178 

0.0757 

Hole 

Size 

20  um 

0.3182 

0.1960 

15  um 

0.8884 

0.6117 

10  um 

4.1877 

2.7428 

7  um 

13.0767 

8.4066 

Fig.  3-34:  Effect  of  the  hole  size  on  the 
damping  ratio.  Decreasing  the  hole  size 
can  be  seen  to  increase  damping. 


Damped  Response 

Here  we  used  the  results  of  the  previous  section  to  model  the  effect  of  the 
damping  on  the  transient  response  of  the  device,  when  subjected  to  a  half-sine  shock 
pulse.  By  including  damping  in  the  model,  we  can  see  that  decreasing  the  size  of  the  hole 
has  a  large  effect  on  the  maximum  deflection  of  the  device. 


Effect  of  hole  size,  damped  response,  0.3  ms,  10g's 


Fig.  3-35:  Effect  of  the  hole  size  on  the  transient  response  of  the  G-Sensor,  when 
subjected  to  a  half  sine  shock  pulse. 

We  are  also  interested  in  determining  the  shock  pulse  which  will  cause  the  G-sensor  to  contact  the 
substrate.  To  this  end,  we  determine  the  shock  level  at  which  the  device  will  touch  the  substrate. 


Table  and  Fig.  show  that  decreasing  the  hole  size  has  a  significant  effect  on  the 
shock  resistance  of  the  G-Sensor.  It  also  shows  that  the  larger  gap  size  will  contact  the 
substrate  at  a  lower 


Gap 

Hole  Size 

2  pm 

4  pm 

25  pm 

23.95  g 

19.83  g 

20  pm 

28.86  g 

25.15  g 

15  pm 

45.82  g 

42.96  g 

10  pm 

137.08  g 

102.80  g 

7  pm 

397.61  g 

265.25  g 

Table  3-15:  Shock  level  that  causes  the  G-sensor  to  contact  the  substrate. 


Shock  threshold 


Fig.  3-36:  Shock  level  that  causes  the  G-sensor  to  contact  the  substrate. 


Recommendations 

Since  it  may  not  be  possible  to  remove  the  holes  completely,  their  size  should  be 
decreased  to  as  small  as  technologically  possible,  with  the  total  number  of  holes 
remaining  the  same.  This  will  allow  the  device  to  be  released  from  the  substrate  and  will 
maximize  damping.  In  addition,  the  increase  in  mass  that  occurs  when  the  hole  area  is 
decreased  will  lead  to  a  more  sensitive  G-sensor.  Adding  mechanical  stops  or  surface 
coatings  to  prevent  stiction  are  other  alternatives;  however  that  is  beyond  the  scope  of 
this  work. 
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Sub  Task  4:  Numerical  study  of  stiction  phenomenon  in  MEMS 

Performer:  Timothy  Singler 
Organization:  SUNY  Binghamton 

Introduction 

Stiction  plays  an  important  rule  in  the  manufacturing  reliability  and  the 
performance  of  MEMS  and  NEMS  devices.  The  small  characteristic  length  scales  of 
micro  devices  emphasize  the  effect  of  surface  forces  over  the  volume  forces  which  in  turn 
results  in  stiction  problems.  Surface  forces  of  various  origins  can  contribute  to  stiction; 
these  include  capillary.  Van  der  Waals,  electrostatic  and  Casmir  forces.  Capillary  forces 
result  from  the  presence  of  water  traces  that  appear  within  the  structure  during 
manufacturing  processes  or  as  a  result  of  capillary  condensation  of  water  from  the  humid 
environment.  The  Van  der  Waals  forces  could  be  attractive  or  repulsive  depending  on  the 
separating  distances  between  different  surfaces  of  the  structure.  When  the  stiction  forces 
become  larger  than  the  actuating  force  for  the  structure,  adhesion  of  the  surfaces  can 
occur  resulting  in  fabrication  failure  or  mat  function  of  the  MEMS  device. 

Recently,  a  great  deal  of  research  has  been  directed  at  understanding  the 
fundamentals  of  stiction  and  at  constructing  models  that  describe  stiction  behavior.  The 
stiction  energy  was  first  measured  by  [1],  by  using  an  array  of  cantilevered  beams  as  a 
test  device.  [2]  studied  the  stiction  of  polysilicon  microbeams  in  a  controlled  humidity 
environment.  They  used  supercritical  CO2  drying  for  freeing  the  beams.  They  found  that 
the  supercritical  CO2  drying  is  able  to  free  a  beam  of  length  up  to  2  mm  length.  Also, 
they  found  that  the  surface  interaction  energy  increases  exponentially  with  the  relative 
humidity.  [3]  derived  a  stiction  model  for  suspending  polysilicon  microbeams  that 
includes  the  effect  of  the  residual  stress  gradients  and  the  post  release  temperature.  They 
used  three  different  materials  for  the  substrate:  bare  silicon,  unmodified  smooth 
polysilicon  and  grain  hole  polysilicon.  They  found  that  the  grain  hole  polysilicon 
substrate  tends  to  decrease  the  effect  of  the  residual  stress  gradient  on  stiction  as  a  form 
of  lower  surface  interaction  energy.  [4]  discussed  the  theoretical  limit  on  the  freestanding 
length  of  cantilevered  beams,  taking  into  account  the  acceleration,  Casmir  and  coulomb 
forces.  They  tabulated  a  set  of  equations  to  determine  the  maximum  length  of  the  beam 


for  which  no  stiction  failure  could  occur.  [5]  studied  the  contact  deformation  of  a  plate¬ 
like  clamped/clamped  micro  beam  and  derived  a  relation  between  the  surface  interaction 
energy  and  the  bending  moment  at  the  edge  of  the  contact  zone,  they  also  found  a  lower 
limit  for  the  gap  between  the  micro  structure  and  the  substrate  and  this  limit  is 
proportional  to  the  structure  length  and  to  the  square  root  of  the  surface  interaction 
energy.  [6]  found  that  for  coated  hydrophobic  beams,  adhesion  is  independent  of  RH  up 
to  a  threshold  value  that  depends  on  the  coating  used.  [7]  and  [8]  introduced  a  physical 
model  for  investigating  the  sensitivity  of  microstructures  to  stiction.  The  model  included 
the  capillary  energy,  Van  Der  Waals  forces,  electrostatic  forces  and  H-bond  forces.  It  also 
took  account  of  the  roughness  of  the  surface  by  assuming  a  Gaussian  distribution  of  the 
separation  distance  between  the  two  surfaces.  [9]  studied  the  effect  of  the  surface  tension 
on  the  performance  of  a  SOI  (silicon  on  insulator)  beam  and  discovered  that  the  adhesion 
forces  are  directly  proportional  to  the  surface  tension  of  the  last  liquid  used  before  drying. 


Background 


Equilibrium  interfaces 


Capillarity 

Capillarity  concerns  the  motions  and  shapes  of  interfaces  in  response  to  the 
presence  of  surface  forces.  An  equilibrium  capillary  surface  refers  to  a  surface  (either 
liquid-liquid  or  liquid-gas)  whose  shape  is  dictated  by  the  equilibrium  between  surface 
and  body  forces.  The  formation  of  equilibrium  capillary  surfaces  requires  an  amount  of 
free  energy  that  is  ultimately  stored  in  the  interface  itself.  Any  change  in  the  equilibrium 
shape  will  result  in  a  change  in  the  equilibrium  free  energy  of  the  interface;  thus  an  equal 
amount  of  work  will  be  needed  to  be  performed  on  the  interface  in  order  to  change  its 
equilibrium  free  energy  by  the  same  amount.  The  equilibrium  shape  of  the  interface  is 
related  to  the  pressure  difference  A P  across  the  interface  by  the  Young-Lap  lace 
equation, 


a p  =  n- 


0) 


where  Rf  and  R2  are  the  principal  radii  of  curvature,  yh.  is  the  surface  tension  of  the 
liquid. 

Effect  of  vapor  pressure  change  on  the  equilibrium  of  interfaces 

A  change  in  the  vapor  pressure  of  the  liquid  will  result  in  a  change  in  the  pressure 
drop  across  the  interface.  This  in  turn  will  result  in  a  change  of  the  equilibrium  shape  of 
the  interface.  The  relation  between  the  vapor  pressure  and  the  equilibrium  shape  of  the 
interface  can  be  derived  from  equilibrium  thermodynamics  and  is  known  as  Kelvin's 
equation: 

R gT\n[— \=Vmyi—  +  — )  (2) 

*  (p°)  U  K2) 

where  Rk  is  the  gas  constant,  P°  is  the  vapor  pressure  of  the  liquid  and  Vm  is  the  molar 

volume  of  the  liquid.  Kelvin's  equation  can  be  used  to  express  the  relationship  between 

P 

the  relative  humidity  RH  and  the  change  in  the  equilibrium  shape.  Since  RH  =  — , 
Kelvin’s  equation  becomes 

«,rln(AH)=K,r,|i-  +  -lj  (3) 

This  may  be  specialized  to  a  spherical  section  of  L/V  interface  by  letting  R  =  R}=  R2 
which  leads  to 

07/ 

RgT\n{RH)==^-  (4) 

A 

Solid  liquid  vapor  interfaces 

Consider  a  static  contact  line  formed  by  a  liquid/vapor/solid  triple  junction  and 
characterized  by  a  contact  angle  9c .  Consider  the  virtual  infinitesimal  displacement  of 

that  contact  line  such  that  liquid  displaces  vapor  and  there  is  an  increase  A A  in  the 
solid/liquid  (S/L)  interfacial  area  and  a  commensurate  decrease  in  the  solid/vapor  (S/V) 
interfacial  area.  The  associated  change  AG  in  the  Gibbs  free  energy  of  the  system  is  then 
given  by 


(5) 


AG  =  A  A  {ys]  -  yJ+AAylv  cos(^) 

Since  the  triple  junction  region  is  in  equilibrium,  AG  and  Young's  equation  follows: 

(r.1  -  r„ )  +  Yh cos(4 )=°  =>r„~ Ysi  =  Yw cos(&c )  ( 6) 

Consider  a  solid  surface  of  fixed  area  As  and  assume  that  an  area  Asl  of  As  is  wetted  by 
a  liquid.  We  emphasize  here  that  Asl  is  capable  of  variation  while  4,  is  a  constant.  The 
total  Gibbs  free  energy  G  associated  with  the  surface  is  given  by 

G  =  (As~  A, k,  +  AiYsi  =  aY„  +  4/ (Ys,  ~  Y sv )  =  A ysr  -  As,yh.  cos (0e )  (7) 
where  use  has  been  made  of  Young's  equation.  Assuming  that  the  S/V  interfacial  energy 
is  constant  so  that  the  term  Aysv  is  also  constant,  the  Gibbs  free  energy  of  the  surface 
may  be  written  as 

G  =  G0-AsiY lycos(0c)  (8) 

This  is  because  the  total  energy  G  can  only  be  defined  to  within  an  arbitrary  constant, 
and  the  constant  term  Aysv  may  be  absorbed  by  that  constant. 

Surface  interaction  energy  measurement 

Mastrengelo  et.  al.[l]  proposed  a  novel  technique  to  measure  the  surface 
interaction  energy.  The  measuring  technique  is  applied  using  an  array  of  cantilevered 
beams  of  different  length  whereby  a  force  is  applied  to  each  beam  until  it  contacts  the 
underlying  substrate.  The  bending  (elastic)  energy  UE  of  an  individual  beam  can  be 
calculated  if  the  unadhered  length  s  of  the  beam  is  known.  Minimizing  the  total  energy 
UT  =  UE  +  T(/  -  5)  with  respect  to  s  then  yields  the  surface  interaction  energy  density 

w  ds 


Analytical  support  for  the  above  measurement  technique  requires  solution  of  the  elastic 
beam  equation 


(9) 


subject  to  the  boundary  conditions 


u(x  =  0)=  0,u'(x  =  0)=  0  ,u(x=  s)  =  -h,n’(x  =  5)=  0  (10) 

here  u  is  the  deflection  of  the  beam  at  the  axial  position  x,  h  is  the  seperation  distance 
between  the  beam  base  and  the  substrate,  and  w  is  the  beam  width.  The  resulting  solution 
for  the  shape  of  the  deformed  beam  is 


u(x)  =  h 


3-2- 


(11) 


The  elastic  energy  may  be  calculated  from  the  above  solution  by  integrating  the 
deflection  equation  as  follows: 

El  N2 


Ue=- 


r^i  dx 

//v2 


dx 


(12) 


where  E  is  modulus  of  elasticity  and  I  is  the  second  moment  of  inertia. 

Another  possible  configuration  for  the  beam  deflection  is  the  arc  shape  (Fig.  1)  in  which 
the  beam  contacts  the  substrate  only  at  its  tip  De  Boer  et  al.[6]  investigated  both 
configurations  and  reached  the  conclusion  that  the  s-shape  beam  is  sufficient  for  the 
study  of  stiction  problem  in  microbeams. 


Surface  Evolver 

The  Surface  Evolver  is  an  interactive  program  for  the  study  of  surfaces  shaped  by 
interfacial  and  other  energies.  A  surface  is  modeled  as  a  union  of  triangular  elements.  The 
program  has  no  user  friendly  interface,  and  an  initial  surface  must  be  defined  in  a  data 
file.  This  surface  is  then  evolved  towards  minimal  energy  by  the  method  of  gradient 
descent. 


Modeling 

Stiction  during  Manufacturing  processes 

The  primary  objective  of  the  current  work  is  to  study  the  effect  of  capillary 
stiction  on  the  operation  of  a  thermally  actuated  switch.  The  switch  is  actuated  by  the 
ohmic  heating  of  a  supporting  beam,  the  resulting  deformation  of  which  actuates  a  central 


beam  element  to  complete  either  an  electrical  or  optical  coupling.  A  micrograph  of  the 
switch  is  shown  in  Fig.  4-1. 


Fig.  4-1.  Micrograph  of  thermally  actuated  switch. 

In  principle,  stiction  could  occur  between  any  and  all  of  the  beam  structures  that 
comprise  the  switch  mechanism.  However,  due  to  its  complexity,  a  Surface  Evolver 
model  of  capillary  stiction  between  all  elements  of  the  complete  switch  and  the 
underlying  substrate  is  formidable.  Instead,  we  focus  on  a  Surface  Evolver  model  of  the 
primary  actuation  element,  the  central  beam.  Two  general  models  of  this  beam  are 
considered:  (1)  a  rigid  beam  with  either  a  round  or  flat  tip  and  which  is  inclined  with 
respect  to  the  substrate;  (2)  an  elastic  beam  which  deforms  due  to  stiction  forces.  In  both 
models,  the  dimensions  of  the  beam  can  be  varied.  Surface  Evolver  models  were  built  for 
four  specific  cases  of  the  two  general  models  described  above.  These  are  described 
below. 

Case  1 

This  case  (Figure  4-2)  is  a  rigid  beam  with  flat  tip.  The  effects  of  angle  of  inclination  and 
beam  tip-substrate  distance  on  the  stiction  force  and  surface  interaction  energy  were 
explored. 


Fig.  4-2.  Case  1 :  Stiction  model  for  an  inclined  rigid  beam  with  straight  tip. 


Case  2 

This  case  (Figure  4-3)  is  a  rigid  beam  with  curved  tip.  The  effects  of  angle  of  inclination 
and  beam  tip-substrate  distance  on  the  stiction  force  and  surface  interaction  energy  were 
explored. 


Fig.  4-3.  Case  2:  Stiction  model  for  an  inclined  rigid  beam  with  curved  tip. 


Case  3 


This  case  (Figure  4-4)  is  a  rigid  beam  with  curved  tip  with  zero  angle  of 
inclination.  Figure  5  shows  a  schematic  drawing  for  this  case.  The  liquid  is  pinned  along 
curve  2  and  along  the  sides  of  the  beam;  the  contact  line  is  free  along  curve  1.  The 
stiction  is  characterized  by  the  parameters  F.  (force  in  z-direction),  F  ,  (force  in  y- 

direction)  and  U A  (surface  interaction  energy);  the  horizontal  force  component  Fx  is  zero 
by  symmetry.  The  force  components  are  calculated  by  the  method  of  virtual 
displacement.  The  surface  interaction  energy  UA  is  defined  as  the  total  interfacial  energy 
formed  between  condensed  phases.  The  liquid/vapor  interfacial  energy  (surface  tension) 
is  denoted  by  yh  and  the  liquid  volume  as  V .  The  contact  angles  of  the  liquid  on  the 
beam  and  substrate  materials  are  denoted  by  0b  and  6S ,  respectively.  The  geometry  of 
the  beam  is  characterized  by  its  width  w  and  its  height  h  above  the  substrate. 


Fig.  4-4.  Case  3:  Stiction  model  for  a  rigid  beam  with  curved  tip. 

There  are  three  relevant  interfacial  areas,  the  solid/liquid  interfacial  area  A*  between  the 
beam  and  liquid,  the  solid/liquid  interfacial  area  Ass,  between  the  substrate  and  liquid  and 
the  liquid/vapor  interfacial  area  A/y.  The  total  interfacial  energy  is  the  sum  of  the 


liquid/vapor  interfacial  energy  and  the  interfacial  energies  of  the  condensed  phases 
(surface  interaction  energies).  Thus 

E-r  =  Eh  +  Esl  +  Esl  (13) 

Making  use  of  equation  (6),  this  expression  may  be  expressed  as 

et  =  l  rA  “  \A,  Yw  cos(^  )dAf,  -  yh  cos(0s  )dAsd  ( 1 4) 

where  ET  is  the  total  energy  of  all  interfacial  regions,  Eh,  is  the  liquid  vapor  interfacial 
energy,  Ed  is  the  interfacial  energy  of  the  beam/liquid  interface  and  Ed  is  the  interfacial 
energy  of  the  substrate/liquid  interface.  Surface  Evolver  provides  a  numerical 
approximation  to  (14).  Since  Surface  Evolver  provides  ET,  the  surface  interaction  energy 
can  be  calculated  as 

E Et  -  Y iv A ,y  ( 1 5) 


Fig.  4-5.  Schematic  drawing  of  case  3 


Figures  4-6  ~  4-8  show  the  variations  of  F_  ,  Fy  and  U  A  for  a  beam  of  width  50 
/urn.  F_  may  be  seen  to  increase  as  V  increases.  This  is  due  to  the  increase  in  the  length 
of  the  contact  lines  along  the  straight  edges  of  the  beam.  F.  also  increases  as  6b 


increases  (with  9S  held  fixed)  in  the  range  0  <9b<  —  .  Fy  is  not  affected  by  the  change 

.  71 

in  V  but  decreases  as  9b  increases  in  the  range  0  <  9b  <  —  .  This  is  because  the  shape  of 

the  curved  interfaces  (curves  1  and  2  in  Fig.  5)  do  not  change  with  V  provided  that  V  is 
sufficiently  large.  For  this  reason,  Fy  is  insensitive  to  changes  in  V.  U  A  increases  as  U  A 

71 

increases  and  decreases  as  9b  increases  in  the  range  0  <9b<  —  \  the  latter  behavior  is 


expected  since  E*  =-  \  Bylvcos(0b)dA* .  The  increase  in  UA  with  V  is  due  to  the 

•*4.3 


increase  in  A*  with  increasing  V. 


Fig.  4-6.  F_  versus  9b  and  V  for  a  beam  of  width  50  fjm . 


Fig.  4-7.  Fy  versus  9b  and  V  for  a  beam  of  width  50  jum. 


Fig.  4-8.  U A  versus  0b  and  V  for  a  beam  of  width  50  fm. 


Figures  4-9  -4-1 1  show  the  variations  of  F:  ,  Fy  and  U A  for  a  beam  of  width  90  fum. 
By  comparing  this  set  of  figures  to  the  previous  set,  it  can  be  seen  that  as  the  width  of  the 


beam  increases  for  a  fixed  value  of  V,  F.  decreases.  This  is  due  to  the  decrease  in  the 
length  of  the  contact  lines  along  the  straight  edges  of  the  beam.  F  increases  as  the  width 

of  the  beam  increases.  This  is  due  to  the  increase  in  the  length  of  the  contact  lines  along 
curves  1  and  2.  U A  increases  as  the  width  of  the  beam  increases  This  is  due  to  the 

increase  in  Af, . 


Fig.  4-9.  F.  versus  0b  and  V  for  a  beam  of  width  90  /jm. 


Fig.  4-10.  F  versus  9b  and  V  for  a  beam  of  width  90  /jm. 


Fig.  4-11.  UA  versus  9b  and  V  for  a  beam  of  width  90  /m. 


Case  4 

Figure  4-12  shows  an  elastic  beam  in  contact  with  the  underlying  substrate  due  to  a 
capillary  bridge.  This  model  has  been  widely  used  to  investigate  the  effect  of  stiction  on 
the  performance  of  MEMS  devices.  Most  of  the  existing  literature  relies  on  crude 
approximations  to  the  surface  interaction  energies.  We  will  exploit  this  model  and  use 
Surface  Evolver  to  accurately  calculate  the  forces  and  surface  interaction  energies.  By 
including  the  elastic  energy  of  the  beam,  we  will  be  able  to  accurately  calculate  the  total 
system  energy.  By  minimizing  this  energy,  we  will  be  able  to  determine  the  length  of 
contact  between  the  beam  and  underlying  substrate.  By  determining  the  critical  system 
parameters  that  insure  that  there  is  no  beam/substrate  contact,  we  hope  to  identify  a  safe 
design  space  to  avoid  stiction. 


Fig.  4-12.  Case  4:  Stiction  effect  on  flexible,  cantilevered  beam. 


In  Use  stiction 

Section  problem  occur  in  MEMS  either  during  Manufacturing  processes  or  during 
the  operation  of  MEMS  devices.  The  former  is  known  as  in-use  stiction.  The  reliability  of 
MEMS  device  can  be  improved  making  use  of  anti-stiction  coating.  Anti-stiction  coating 
usually  applied  to  MEMS  surfaces  in  order  to  improve  the  surface  characteristic  like 
surface  energy,  hydrophobiciy  and  surface  roughness.  The  successful  coating  can  achieve 
by  improving  the  characteristics  in  such  away  the  surface  will  have  high  contact  angle, 
low  Van  der  Waals  interaction  energy  and  acceptably  high  surface  roughness.  In  the 
current  work,  a  mathematical  model  has  been  built  in  order  to  achieve  parametric  studies 
of  in-use  stiction  phenomena.  In-use  stiction  occur  when  the  adhesion  force  is  greater 
than  the  restoring  force  of  the  structure  or  when  the  adhesion  energy  is  considerably  high. 
In  the  current  model,  stiction  phenomena  is  characterized  by  the  adhesion  force  and  the 
adhesion  energy.  The  effect  of  different  parameters  on  stiction  likelihood  has  been 
studied.  The  parameters  are  relative  humidity,  contact  angles,  surface  roughness, 
Hamaker  constant  (associated  with  Van  der  Waals  interaction),  modulus  of  elasticity  and 
Hardness  of  MEMS  materials. 


Mathematical  formulation 


Surface  roughness 


Fig.  4-13.  Rough  and  smooth  surfaces  in  contact 

As  shown  in  figure  4-13,  Two  rough  surfaces  in  contact  can  be  modeled  as  one  surface  is 
rough  and  elastic  and  the  other  surface  is  rigid  and  smooth  [10].  The  surface  roughness  of 
the  surface  is  modeled  following  Greenwood  and  Williamson  model  [11],  the  surface 
roughness  has  been  modeled  assuming  that  all  the  asperities  have  the  same  radius  of 
curvature  and  the  surface’s  heights  distribution  is  Gaussian.  The  properties  of  the 
equivalent  rough  surface  are  as  follows 

a2  =  of  +  o\  (16) 


Where: 

a  is  the  standard  deviation  of  the  equivalent  rough  surface  height  distribution, 
er,  is  the  standard  deviation  of  the  first  surface  height  distribution. 

<r2  is  the  standard  deviation  of  the  second  surface  height  distribution. 

E'  is  the  equivalent  surface  modulus  of  elasticity. 

£■,  is  the  first  surface  modulus  of  elasticity. 

E2  is  the  second  surface  modulus  of  elasticity, 
uis  Poisson  ratio. 


Whitehouse  et.  al.  [12]  proposed  a  probability  distribution  function  for  asperity’s  peak 
height  distribution  and  they  introduce  a  formula  for  the  statistical  property  of  the  asperity 
peak  distribution. 


1  +  erf\ 


2 


(18) 


Where  p  is  function  of  the  sampling  interval  (1)  and  the  spatial  autocorrelation  /?’ 


(19) 


op  =  0.9<t 

(20) 

20cr 

(21) 

OpRrj  -  constant 

(22) 

Where:  op  is  the  standard  deviation  of  the  peak  height  distribution,  R  is  the  mean 
peak  radius,  77  is  the  peak  density. 

McCool  et.  al.  [13]  proposed  a  relation  between  a ,  ob,  R  and  77 


a 


2 


=  op  + 


3.71 7eT4 
77  2R2 


(23) 


With  cr  =  0.9 a ,  one  can  solve  for  abRrj . 


obRtj  =  40.2e  3 


(24) 


Structural  forces 

When  two  surfaces  have  contact  their  asperities  will  indent  each  other  which  will  result  in 
a  reaction  forces  that  try  to  resist  the  indentation  process.  Since  the  asperities  peak  is 
assumed  to  have  a  constant  radius  of  curvature  then  Hertzian  contact  theory  can  be 
applied. 


Fig.  4-14.  Elastic  sphere  indenting  smooth  and  rigid  surface. 

The  literature  contains  empirical  formulae  for  the  elastic,  elasto-plastic  and  fully  elastic 
forces  of  the  Hertzian  contact  [14-15].  The  structural  force  P(z)  is 
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is  a  function  of  the  Poison  ratio  and 


(26) 


Where  H  is  the  hardness  of  the  surface  and  K'  is  a  constant  depends  on  the  Poisson  ratio 
K  =  0.454  +  0.4 \v. 


The  total  reaction  force  per  nominal  surface  area  applied  on  the  rough  surface  is 
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Capillary  Condensation 
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Fig.  4-15.  Liquid  meniscus  between  sphere  and  half  space. 

Orr  at.  al.  [16]  studied  a  liquid  meniscus  between  sphere  and  half  space  where  the 
sphere  is  just  touching  the  half  space  without  indentation.  In  the  current  analysis  Orr.  at. 
al.  analysis  has  been  generalized  to  include  all  possible  configurations  of  the  separating 
distance  between  the  sphere  and  the  half  space. 


Figure  4-15  shows  the  geometry  of  sphere  and  half  space  that  they  have  a 
separation  distance  D  between  them.  Where  D  =  z  -z ,  rl  is  the  meridional  radius,  r2  is 

the  azimuthal  radius,  <9,  is  the  contact  angle  between  the  liquid  and  the  sphere’s  surface, 
92  is  the  contact  angle  between  the  liquid  and  the  flat  surface,  cp  is  the  filling  angle  of  the 
liquid  over  the  sphere’s  surface,  R  is  the  radius  of  the  sphere  and  zm  is  the  maximum 
height  of  the  liquid  meniscus. 


By  analyzing  the  geometry  of  Fig.  4-15,  ones  can  find  that  The  azimuthal  and  the 
meridional  radii 

9?(l  -  cosip)  +  D 


cos(9,  +^)+cos(#2) 
7?sin(#>) 
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2  sin(6l1  +  <p ) 

Given  the  contact  angles,  the  separation  distance  D  and  the  radius  of  the  sphere  R, 
ones  can  solve  for  the  filling  angle  making  use  of  the  Kelvin’s  equation  [17-18] 

RT 
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Where  Rg  is  the  real  gas  constant,  T  is  the  temperature  of  the  air,  V,  is  the  liquid 
molar  volume,  RH  is  the  relative  humidity  of  the  air  and  y  is  the  surface  tension  of  the 
liquid. 

The  meniscus  forces  are  consist  of  the  surface  tension  forces  and  the  Lap  lace 
pressure  force. 
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eq 


The  adhesion  energy  is  defined  as  the  difference  between  the  interfacial  energies  of  the 
state  of  contact  and  the  state  of  infinite  separation  [19] 

£,  =  r„iA„,  +  r^A,2  +  yK  ~rcos{9])Asll  -ycos{92)Asl2 

=  YsrArl  +  Ysv2^sv2 


(32) 

(33) 


ec(z)  =  E2-Et  (34) 

ec(z )  =  ycos  (0,)As/l  +  /cos  (02)Asl2-yAly  (35) 

Where  ysv]  is  the  interfacial  energy  of  solid  vapor  interface  of  the  sphere  surface, 
y^,2  is  the  interfacial  energy  of  the  solid  vapor  interface  of  the  substrate,  Asv]  is  the  total 
surface  area  of  the  sphere,  Asr2  is  the  total  surface  area  of  the  substrate,  ^s/)is  the 
solid/liquid  interface  between  the  liquid  and  the  sphere  surface,  Asl2  is  the  solid/liquid 
interface  between  the  liquid  and  the  substrate  surface  and  Ah,  is  the  liquid/vapor 


interfacial  energy. 
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Ah  =  2 n  Jr, (sin(^)  +  sin(6»,  +  <p)) -  ^r,2  -  (z  -  zeq  +  r,  cos(6'1 )) dz  (38) 


A,v  =  2  nr? 


(cos(6)1 )  +  cos(#,  +  £>)Xsin(#>)  +  sin(#,  +  $?))  +  ^  +  ^  +  9 
--^-sin(2<9,  +2^)-isin(26>2) 


The  Capillary  force  and  adhesion  energy  is  follows 


co 

FcoP(Zeg)=Tl  jfJZMZ)^ 

—  00 

m 

KpkJ=r7  feAz)0(z)dz 


F'op  and  E'cap  have  values  different  from  zero  only  when  the  meniscus  is  stable. 
The  meniscus  is  considered  instable  when 

r,sin(^,  +  (p)+  7?sin(<j?)<  r,  (42) 


Van  der  Waals  interaction 


The  derivation  of  Van  der  Waals  interaction  energy  of  sphere  above  a  half  space 
without  the  existence  of  liquid  meniscus  can  be  found  in  [18].  In  the  current  analysis. 
Van  der  Waals  interaction  energies  and  forces  have  been  derived  for  the  cases  of  sphere 
indenting  half  space  without  the  existence  of  liquid  meniscus,  A  sphere  above  half  space 
with  existence  of  liquid  meniscus  and  a  sphere  indenting  a  half  space  with  the  existence 
of  liquid  meniscus. 


(4,-4)* 

2 yc+Di  4/U+4 

Jc+D  (Ao-A/)D_ 

(4.-4,)* 

~2 X+D^r+Iftd-I) 

Jc+D  (4-4/K 

4/ 

6D 

A£ 

M0 

Z^Zeg-do’<P^0 

Z^Zeg-do^0 

z<zev-do,<p=0 

Z^Zeg-do’<P=0 


(43) 


f(z)= 


XAa-A)R 

\+D{  Ak+rf' 

Jr+D  (Ao~A,)D_ 

( Aa~Al)R 

jc+D  {Ao-aK 

Al 

6/T 

AJR 

Ml 

z<zeq-do,(p>0 

Z>Zeg~do’<P^ 

Z^'g-Clo^O 

Z^Zeg-doW=® 


(44) 


Where 

yc  =  r,(cos(0,  +  <p)  +  cos(#2))-  D 

(45) 

D  =  z  -D 

eg 

(46) 

do  is  the  cutoff  distance  of  the  attractive  Van  der  Waals  interaction  d0=.  2  nm 

Aha  is  the  Hamaker  media  of  the  materials  of  the  sphere  and  the  half  space  interaction  in 

air. 

Ahw  is  the  Hamaker  media  of  the  materials  of  the  sphere  and  the  half  space  interaction  in 


water. 


Hamaker  constant  can  be  found  using  Lifshitz  theory  and  the  combining  relation  [19] 


-2m 

4 


,*1+^3  J 


+ 


3 hve  (»,2-«32)? 

16V2  („*  +  nf¥2 


(47) 


An]  is  Hamaker  constant  of  two  media  of  material  1  interacting  in  medium  of  material  3, 
s,  is  the  dielectric  constant  of  material  1,  e3  is  the  dielectric  constant  of  the  separating 


media,  «,  is  the  index  of  refraction  of  material  1,  n3  is  the  index  of  refraction  of  material 
3,  k  is  the  Boltzmann’s  constant,  T  is  the  temperature,  h  is  Blanck’s  constant,  ve  is  the 

Plasma  frequency  typically  equal  to  3e15Hz  [4].  The  values  of  dielectric  constant  and 
index  of  refraction  are  shown  in  table  4-1 . 


Material 

£ 

n 

air 

1 

1 

Water  [4] 

80 

1.333 

Silicon 

11.8 

HF 

83.6 

1.1574 

Si3NA 

4.2 

Si02 

4.41 

IBH 

Table  4-1 


Table  4-2  shows  values  of  Hamaker  constant  for  different  material.  Hamaker 
constants  in  table  2  can  be  found  in  the  literature  as  indicated  by  the  reference  number. 
Also,  Equation  47  is  used  to  calculate  Hamaker  constant  when  no  reference  is  cited. 
Hamaker  constants  in  table  3  are  found  by  considering  that  the  energy  of  adhesion  is 
determined  by  the  properties  of  the  adsorbed  film  even  when  it  is  only  monolayer 
thick[19]. 


A 

n SiHFSi 

625  zJ 

A 

■ri5i02  nater  Siy  N4 

19  zJ  [20] 

A 

^HF  water  HF 

15.9  zJ 

A FDTSairFDTS 

50  zJ  [21] 

A 

'nSiSi02  Si 

A 

^FDTS  water  FDTS 

<4.6  zJ  [22] 

A 

n Si02  air  Si A 

103.8  zJ  [20] 

A 

^SiOx  water  Si O 2 

4.3  zJ  [23] 

A 

**octadeceneair  octadeeene 

50  zJ  [19] 

A 

** octadeeene  water  octadeeene 

5  zJ  [19] 

A 

yi OTSairOTS 

190  zJ  [24] 

A 

Si  air  Si 

690.71  zJ 

A 

voter 

37  zJ  [19] 

A SiO:  otr  S,0  , 

64.4  zJ 

Tabl 

e  4-2 

Numerical  Solution 

There  are  two  solution  domain,  the  first  domain  is  the  heights  of  the  asperities 
peak  z  and  the  second  is  the  location  of  the  equilibrium  plane  zeg.  To  achieve  grid 

independent  solution.  The  node  numbers  on  the  z-domain  is  selected  to  be  2000  nodes 
over  the  range  (-  R  <  z  <  10/?).  The  applied  forces  applied  on  the  rough  surface  are 
calculated  as  in  Eqns.  (1 1,25,31).  The  numerical  solution  is  implanted  by  solving  for  the 
static  equilibrium  using  secant  method  with  a  convergence  criterion  of  £<10  8 .  The 
solution  procedure  is  required  two  initial  guesses  of  z  in  order  to  start  the  process  of  the 
root  finding  according  the  following  equation. 

Where: 


zi+'=z'  -- 

eg  eg 


(47) 


G(0=  Fe(ZJ~  Kp(Z  J-  FL(ZJ 


/+1  i 

z  z 
£  e<! 

4 


(48) 

(49) 


Kelvin’s  equation  Eqn.  (30)  is  needed  to  be  solved  in  order  to  find  filling  angle  of 
the  liquid  meniscus  and  the  capillary  forces.  Since  Eqn.  (30)  has  multiple  roots,  it  was 
solved  using  bisection  method.  The  interesting  solution  domain  of  Eqn.  (4)  is 


0<<p<y2.  After  solving  Kelvin's  equation,  the  applied  forces  are  calculated  and  then  a 

numerical  solution  of  the  static  equilibrium  plane  is  found.  After  that,  the  energy  of 
adhesion  is  found. 


E m=E  +E 

adh  cap  van 


(50) 


Results  and  Discussion 

The  literatures  contain  a  plenty  of  papers  that  study  the  energy  of  adhesion 
experimentally,  i.e.  (21-24).  The  experimental  method  of  measuring  the  energy  of 
adhesion  was  developed  by  [26-27].  Table  4-3  shows  a  comparison  between  the  energy 
of  adhesion  as  measured  experimentally  in  the  literature,  dry  stiction  model  of  Harriri  et. 
al.  [29]  and  the  value  of  the  energy  of  adhesion  from  the  current  model.  As  a  validation 
for  the  code  the  case  of  oxide  coating  polysilicon  with  zero  contact  angles  and  95% 
relative  humidity  has  been  done.  In  this  case,  the  current  model  estimated  the  adhesion 
energy  to  be  79.4  mJ/m2.  The  adhesion  energy  of  two  smooth  surface  which  have  zero 
contact  angles  is  2 y  =  150.4  mJ/m2. 


Coating 

(deg) 

(deg) 

a, 

nm 

nm 

cr 

nm 

A. 

zJ 

A/ 

zJ 

^adh 

Exp 

pJIm2 

Eadh 

Ref. 

29 

pj/m2 

E adh 

theo 

pj  /  m 

OTS  [22] 

112 

112 

12 

12 

16.97 

190 

59.3 

3 

.068 

.08 

OTS[22] 

112 

112 

3 

3 

4.24 

190 

59.3 

30 

5.44 

6.59 

FDTS  [22] 

115 

115 

12 

12 

50 

ESI 

2.1 

.46 

.014 

FDTS  [22] 

115 

115 

3 

3 

50 

m 

8 

6.44 

Octadecene[24] 

104 

104 

4 

1.89 

| 

Oxide  [23] 

0 

0 

2* 

2 

2.8 

64.43 

4.3 

8000 

562.3 

*  Ref.  [23]  do  not  contain  the  value  of  the  surface  roughness  of  the  oxide  layer  the  value  taken  from  ref.  [25]. 

a  —  A  A  —  A 

** ha  /*coalinfi  air  coaling  '‘hi  '‘couliny  wntcrcoatinp 

Table  4-3 


The  reason  why  the  current  model  adhesion  energy  is  less  than  the  theoretical  one 
is  due  to  the  surface  rough.  The  liquid  can’t  cover  all  the  surfaces  because  of  the  different 
asperities  heights  which  will  result  in  the  existence  of  liquid/vapor  interfaces.  The 
existence  of  liquid/vapor  interfaces  decreases  the  adhesion  energy  as  can  be  predicted 
from  Eqn  (35). 

Figures  4-16  and  4-17  show  the  adhesion  energy  and  the  adhesion  force  as  a 
function  of  the  relative  humidity  and  The  Hamaker  constant  across  the  air.  Other 
parameters  have  the  following  values  Ahl  =  50  zJ,  er=  5  nm,  /?*=500  nm,  <9,  =  #,=30°, 

V,  =  18.86^ m^mol,y,v  =  75. 2<?'3%,  T  =  300  K ,£*  =  163  GPa,  and  H  =  10  GPa.  As 

seen  in  Figs(  16-1 7),  The  Adhesion  energy  and  force  decrease  as  the  relative  humidity 
increases  in  the  range  of  low  relative  humidity.  Since  Hamaker  constant  across  liquid  is 
less  than  Hamaker  constant  across  air,  the  formation  of  small  volume  liquid  menisci  will 
result  in  decreasing  the  adhesion  energy  and  the  adhesion  force.  At  high  relative  humidity 
range,  Capillary  effect  dominates  other  effects.  Therefore,  the  adhesion  energy  and 
adhesion  force  increase  as  the  relative  humidity  increases  in  high  humidity  range. 


Fig.  4-16.  Adhesion  energy  versus  relative  humidity  for  different  values  of  Hamaker 

constant  across  air 


Fig.  4-17.  Adhesion  force  versus  relative  humidity  for  different  values  Hamaker  constant 

across  air 

Figures  4-18  and  4-19  show  the  adhesion  energy  and  the  adhesion  force  as  a 
function  of  the  surface  roughness  and  The  modulus  of  elasticity.  Other  parameters  have 
the  following  values  ,4^=500  zJ ,Ahl  =  50  zJ,  RH  =  70%,  /?*=500  nm,  0,  -d2  =30°, 

V/  =18.8 6e~k  m  /mol  ,yh.  =  75.2e~3  ,  T  =  300  K  and  H  =  10  GPa.  As  seen  in  Figs(  1 8- 

19),  The  adhesion  energy  and  the  adhesion  force  decrease  as  the  surface  roughness 
increases.  As  the  surface  roughness  increases  the  z-range  of  the  surface  heights  will 
increase  and  so  the  liquid  can’t  perfectly  wet  the  rough  surface  even  for  the  case  of  zero 
contact  angle  which  tends  to  decrease  the  adhesion  energy  and  force  of  the  capillary 
action.  Also,  the  wide  z-range  will  result  in  decreasing  the  effect  of  Van  der  Waals 
interaction.  Figs(  18-19)  also  show  the  effect  of  the  modulus  of  elasticity  on  adhesion.  As 
the  modulus  of  elasticity  increases  the  adhesion  energy  and  the  adhesion  force  decrease. 
The  decrease  in  the  adhesion  energy  and  force  with  increasing  modulus  of  elasticity  is 
due  to  the  increase  in  the  rigidity  of  the  surface.  As  the  rigidity  increases,  the  equilibrium 
distance  between  the  two  surfaces  increases  due  to  the  increase  in  the  structural  forces. 
Therefore,  the  adhesion  energy  and  force  decrease  as  the  modulus  of  elasticity  increases. 


Fig.  4-18.  Adhesion  energy  versus  surface  roughness  for  different  values  of  modulus  of 

elasticity 


Fig.  4-19.  Adhesion  force  versus  surface  roughness  for  different  values  of  modulus  of 


Figures  4-20  and  4-2]  show  the  adhesion  energy  and  the  adhesion  force  as  a 
function  of  the  surface  roughness  and  The  hardness.  Other  parameters  have  the  following 
values  ^=500  zl,Ahl  =  50  zJ,  RH  =  70%,  /?*=500  nm,  6,)  =  02=3O°, 

V,  =  \8Me-bm* /mol, ylv  =  75. 2e^N/m,  7  =  300  K  and£*=  163  GPa.  As  seen  in  Figs 

(4-20-4-21),  as  the  hardness  increases  the  adhesion  energy  and  the  adhesion  force 
increase.  This  behavior  is  unpredictable,  since  as  the  hardness  increases  the  structural 
forces  increase  which  will  result  in  decreasing  the  adhesion  energy  and  the  adhesion 
force.  Here,  we  report  this  result  as  it  is  since  we  can't  know  what  causes  the 
unpredictable  behavior  of  the  adhesion  energy  and  force  as  function  of  the  hardness. 


Fig.  4-20.  Adhesion  energy  versus  surface  roughness  for  different  values  of  the  hardness 


Fig.  4-21.  Adhesion  force  versus  surface  roughness  for  different  values  of  the  hardness 


Figures  4-20  and  4-21  show  the  adhesion  energy  and  the  adhesion  force  as  a 
function  of  the  contact  angles.  Other  parameters  have  the  following  values  Aha  =500 


zJ  ,Ah,  =  50  zJ,  RH  =  70%,  /?*  =500  nm. 


a  = 


nm. 


V!  =  1 8.86e-6 mi3 /mol , yh,  =  75.2e~3N/m’  T  =  300  K>£’=  163  GPa  and  H  =  10  GPa.  As 

seen  in  Figs(  4-20-4-21),  as  the  contact  angles  decrease  the  adhesion  energy  and  force 
increase  due  to  the  increase  in  the  wetted  area  of  the  surfaces.  The  adhesion  energy  and 
force  decrease  toward  an  asymptotic  value  as  the  contact  angles  increase.  As  the  contact 

71 

angles  become  larger  than  “the  surfaces  become  hydrophobic  which  will  result  in  an 

imperfect  wetting  of  the  surfaces.  Also,  as  the  contact  angles  the  liquid  menisci  become 
instable.  Consequently,  At  high  contact  angles.  Only  Van  der  Waals  interaction  energy 
and  force  are  exist.  Therefore,  as  the  contact  angles  increase  the  energy  of  adhesion  and 
force  will  decrease  asymptotically  toward  the  values  of  Van  der  Waals  interaction  energy 
and  force. 


Conclusions 

1 .  G-W  surface  roughness  model  is  an  efficient  model  that  can  be  used  to  model  the 
stiction  phenomenon. 

2.  The  advantages  of  G-W  model  is  that  it  gives  the  possibility  for  a  realistic 
modeling  of  the  capillary  condensation  and  its  interaction  energies  and  forces. 
Also,  it  makes  it  possible  to  study  the  effect  of  the  liquid  meniscus  on  Van  der 
Waals  interaction  energies. 

3.  The  Van  der  Waals  interaction  is  the  dominant  cause  of  stiction  at  low  relative 
humidity  while  capillary  condensation  is  the  dominant  cause  of  stiction  at  high 
relative  humidity. 

4.  At  low  relative  humidity  the  adhesion  energy  and  force  decrease  as  the  relative 
humidity  increase  due  to  the  formation  of  liquid  menisci. 

5.  As  the  Surface  roughness  increases  the  adhesion  energy  and  the  adhesion  force 
decrease. 

6.  As  the  contact  angles  increase  the  adhesion  energy  and  the  adhesion  force 
decrease. 

7.  As  the  modulus  of  elasticity  increases  the  adhesion  energy  and  the  adhesion  force 
decrease. 

8. 
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Appendix 

Matlab  scripts  that  are  used  to  study  the  in-use  stiction 

function 

[  FC1 ,  FE1 ,  FV1 ,  EE1 ,  EC1 ,  EV1 ) =rough (gama, thl,th2,Ahl,Ahg,est,H,segl,betal,R 
H,  T,  VI , n ) 

%  A  program  for  estimating  the  stiction  force  taking  into  account  the 
%  capiilary  condensation 

0.5,0.0.0000.0.0.000005.0.0.0,0.5.00.0.0.5.00.0 

O'BOOOOOOOOOOOOOOOaC’O'OOOOOOGO'O 

%  parameters  initiation 

0.0.0.0.00000.0.00000.9,0.9,00.0.0000.0.5,9.0 

O'OOOOOOOOOOOOCOOOOO'OOOOOCOOCC 

%Paramters  difinition 
betal=betal/ segl ; 

R=betal/'2/20; 
est=est*segl/ gama ; 
do=- ,17e-9/segl; 


H  =  H/gama*segl;  %  hardness 

K  =  .545+ .41* .25;  %  related  to  poisson  ration  v=.25 
KK=4/3*est/  (1-  .25A2)  /pi/'2; 

EE=est/ (l-.25^2) /2; 
zhc= (pi*K*H/2/EE) *2*R; 
seg= . 9; 

Ahg=Ahg/gama/segl/'2  ; 

Ahl=Ahl/gama/segl^2; 

zm=.47;  %  mean  of  beak  distribution 

%  R=beta;  %  mean  radius  of  the  curvature  of  the  surface 
eigta=41 . 26e-3/R/seg; 

Rg=8.314;  %  Ideal  gas  constant 
rho=exp(-0.23) 

%%%  descritizing  the  z-dimension 

%  the  domain  between  zhc  and  do 

dzl=  ( 0 . l*R+do) /n; 

dz2=  (110*zhc-do) /n; 

dz3= (R- 110* zhc) /n; 

z(l)=-.l*R; 

for  i=2:n 

z  (i ) =z ( i  —  1 ) +dzl ; 

end 

z (n+1 ) =do; 
for  i=n+2:2*n 

z (i)=z ( i — 1 ) +dz2 ; 

end 

z(2*n+l) =110*  zhc; 
for  i=2*n+2:3*n 

z (i)=z ( i— 1 ) +dz3; 

end 

%  break 

y=-5*R: ( 15*R+R) /10/n ; 15*R; 

%  break 

%  z=-R:dz: 10*R;  %  descritization  in  the  z  direction 

x=z; 

dz=0; 

zeql=0; 

zeq0=0; 

%  zeq0=200e-9; 

%  zeql=zeqO+dz; 

%  err=100; 

%  kkk=l ; 

%  k=0 ; 

%%%%%%%%% 


5-e-9-e-2-9-9:5-0:g;9-9-5.9.9.9-9'5-9-2-5-9:9-9.5-9.9-9:9'9-e-S-^&9- 
oo^oo^^o^^ooooooooooo^ooo^o^  ooo  ooo  o 

%%%  initiation 
W ( 1 : length (x) ) =0; 
fh ( 1 : length (x) ) =0; 
e (1 : length (x) )=0; 
f  (l:length(x) )=0; 
ps2  ( 1 : length (x) ) =0; 

Fc  ( 1 : length ( z ) ) =0; 

Ec ( 1 : length ( z ) ) =0; 

Fe ( 1 : length ( z ) ) =0; 


f cO ( 1 : length ( z ) ) =0; 
f cl ( 1 : length ( z) )  =  0; 


)%%%%%%%%%%%%%%%%%%%%% 


%  building  the  probability  distribution  function 
for  m=l : length ( z) 

phi 11 (m)  =  ( 1/sqrt (2*pi ) *exp ( -1  * ( z (m) -zm)A2/2/segA2) )/seg; 

%  phi 11 (m)  =  (1/sqrt (2*pi ) *exp (-1  * ( z (m)  - 

zm) A2/2/segA2) )/seg*(l+erf((z(m)-zm) /sqrt (2)/seg) )  A2 /4; 

phi 11 (m) = (1/sqrt (2*pi ) *exp (-1* (z(m) ) A2/2) )  * ( 1+erf ( ( z  (m) *sqrt ( ( 1- 

rho) / (1+rho) ) ) /sqrt (2) ) ) A2/4; 

end 

%  break 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%%%  Solving  for  the  filling  angle 

%  ps0=fcapz (x,zeqO,R,thl,th2,T,RH,Rg,Vl/ gama, segl ) ; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


11%%%%%  finding  the  elastic  force 

%  [fe0,ee0]  =  Felastic ( z,  zeqO, R, est , segl , gama) ; 

[ fe, ee, f eu, zhc]  =  Felastic ( z, zeql , R, est , H, segl , gama) ; 
psl=f capz (x, zeql, R, thl , th2 ,  T,  RH, Rg, VI, gama, segl, zhc, dz ) ; 

%%% 

%%%%%%  finding  the  capilary  force 

[ f c, ec) =f capx (x, psl ,  z,  zeql,  R,  gama ,  thl, th2, RH, VI, Rg, T, segl , zhc, dz ) 
%%%%%%%%%%finding  Van  der  Waal  forces 

( fv, ev] =fvdw (Ahl, Ahg, z, zeql,  R,  x, psl , thl , th2, segl, gama, zhc, dz) ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


zeq0= .01; 
kkk=l; 
f el=f e; 
f cl=f c; 
fvl=fv; 
evl=ev; 
ecl=ec; 
eel=ee; 

%  break 

for  m=l :length (z) 

phi  11  (m)  =  (l/sqrt(2*pi)*exp(-l*(z(m)-zm)/'2/2/segA2))/seg; 

%  phi 11 (m)  =  (1/sqrt (2*pi ) *exp (-1* ( z  (m)  - 

zm)A2/2/segA2) )/seg*(l+erf( (z(m)-zm)/sqrt(2)/seg) )A2/4; 

phi 11 (m) = (1/sqrt (2*pi ) *exp (-1* (z(m) )A2/2) )* (l+erf( (z(m) *sqrt ( ( 1 — 

rho) / (1+rho) ))/sqrt(2)))A2/4; 

end 

err=100; 

%  for  k=2 rlength (y) 
k=0  ; 

while  err  >  le-10 
k=k+l ; 
zl=z+zeql ; 
for  m=l : length ( z) 

phil (m) = (1/sqrt (2*pi) *exp (-1* (z(m)-zm)A2/2/segA2))/seg; 

%  phill (m)  =  ( 1/sqrt (2*pi ) *exp (-1* ( z  (m)  - 

zm)A2/2/segA2))/seg*(l+erf((z(m)-zm)/sqrt(2)/seg))/'2/4; 

phil (m)  =  (1/sqrt ( 2  *  pi ) * exp (-1  * (zl (m) )A2/2) )* (l  +  erf( (zl (m)*sqrt ( (1- 

rho) / (1+rho) ) ) /sqrt  (2) ) ) A2/4 ; 

end 

FCl=integ (phil .*fcl*eigta,zl) ; 

FEl=integ (phil .*fel*eigta,zl); 


FVl=iriteg  (phil  .*fvl*eigta,  zl )  ; 


)%% 


"c  "c  "6  "6  "o 


uOCO.OOCCCO.GO.O.t^O.C.C 

3COCOOCOOOXC^>OO^io' 


zl=z+zeqO; 

for  m=l : length ( z ) 

phil  (m)  =  (l/sqrt(2*pi)  *exp  (-1*  (z (m) -zm) "2/2/segA2)  ) /seg; 

%  phi 11 (m)=(l/sqrt(2*pi) *exp (-1* ( z (m) — 

zm)  A2/2/segA2)  )  /seg*  (1  +  erf  (  (z  (m)-zm)/  sqrt  (2)/seg)  )/'2/4; 

phil  (m)  =  ( 1/ sqrt  (2*  pi)  *exp(-l*  (zl (m) ) A2/2) ) *  (1+erf  (  (zl  (m)  *sqrt  (  ( 1  - 

rho) / (1+rho) ) ) /sqrt (2) ) ) ^2/4; 

end 

FC0=integ (phil ,*fcl*eigta, zl) ; 

FE0=integ (phil . *fel*eigta, zl ) ; 

FV0=integ (phil . *fvl*eigta, zl ) ; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


G=integ (phil . *ecl , zl ) *eigta+integ (phil . *evl ,  zl) *eigta+integ( phil.* eel, 
1 ) *eigta; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%? 

FO=FEO-FVO-FCO ; 

F1=FE1-FV1-FC1 ; 


%  %  %  % 


FFO (k) =FE1 ; 

FF2 ( k) =FC1 ; 

FF1 (k) =F1; 

FF3 (k) =FV1 ; 

EE1 ( k) =integ (phil . *ecl , zl ) *eigta; 
EE2 (k) =integ (phil . *evl , zl ) *eigta; 
EE3 ( k) =integ (phi 1 . *eel , zl ) *eigta; 
GG ( k ) =G ; 
zz (kkk) =zeqO; 
kkk=kkk+l ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

if  FI  ~=  FO 

zeq=zeql-Fl* ( zeql-zeqO) / (F1-F0)  ; 
else 


zeq=.75*zeql 

end 

err=abs ( ( zeq-zeql ) /zeq) *100 

zeq0=zeql ; 

zeql=zeq; 

%  zeql=zeq; 


end 

%  inn=find(FFl  <  0) ; 

%  zeql=y (inn (length (1 ))) ; 

%  k=inn (length ( 1 ) )  ; 

zl=z+zeql ; 

for  m=l : length ( z ) 

phil  (m)  =  ( 1/sqrt  (2*pi)  *exp  (-1*  (z(m)-zm)''2/2/segA2)  )/seg; 

%  phill(m)=(l/sqrt(2*pi)*exp(-l*(z(m)- 

zm),v2/2/seg/'2))/seg*(l+erf((z(m)-zin)/sqrt(2)/seg))^2/4; 
phil(m)=(l/sqrt(2*pi) *exp (-1* (zl(m) )/'2/2) )* (l+erf( (zl (m) *sqrt ( (1- 
rho)  /  (1  +  rho)  )) /sqrt  (2)  )  )/'2/4; 
end 

FCl=integ (phi l.*fcl,zl)*eigta; 

FEl=integ(phil.*fel,zl)*eigta; 

FVl=integ (phi l.*fvl,zl)*eigta; 

G=integ(phil.*ecl,zl)*eigta+integ(phil.*evl,zl)*eigta; 


ECl=-l*integ(phil.*ecl,zl)*eigta; 
EVl=-l*integ(phil. *evl , zl ) *  eigta; 
EEl=integ (phil .*eel,zl)*eigta; 


function  ps2=f capz (x,zeq,R,thl, th2, T,  RH, Rg, VI, gama, seg, zhc, dz) 
ps2  ( 1 : length (x) ) =0; 

%%%%%%%%%%%%%%%%%%%%%%%%%% 

xc=0; 

for  i=l : length (x) 
err=100; 
h=zeq-x (i)  ; 
if  h  <=  -R 
h=-R; 

end 

if  h  <=  0 

xa=acos ( 1+h/R) ; 

else 

xa=0; 

end 

xb=pi/2+ . 01 ; 

EPSILON  =  le-8; 
phii0=xa; 

fa  =  (( -1* (cos (th2 ) +cos (thl+phiiO) )/(( 1- 
co s(phiiO) )+h/R)+l/sin (phiiO) ) ) -seg*R*Rg*T*log (RH) /gama /VI; 
phii0=xb; 

fb  =  ( (-1* (cos (th2)+cos (thl+phiiO) ) /( (1- 
cos (phiiO) ) +h/R) +l/sin (phiiO) ) ) -seg*R*Rg*T*log (RH) /gama /VI ; 
iter=0; 

while  (  err  >  EPSILON  ) 
iter=iter+l; 

xc  =  abs (  xa  +  xb  )  /  2; 
phii0=xc; 

fc  =  ( (-1* (cos (th2) +cos (thl+phiiO) ) / ( (1- 
cos (phiiO) )+h/R)+l/sin (phiiO) ) ) -seg*R*Rg*T*log (RH) /gama /VI; 


xc=0; 


err=  abs  (  fc  )  ; 
if  (  fc  ==  eps  ) 
xa  =  xc; 

elseif  (  sign(fb) 
xa  =  xc; 
fa  =  fc; 

else 

xb  =  xc; 
fb  =  fc; 

end 

if  iter  >  200 
err=eps; 
if  h  >  0 
xc=0; 

else 

end 

end 


end 

if  xc  >  pi/2 


sign(fc)  <=  0  ) 


xc=pi/2 ; 

end 

ps2 (i)=xc; 

end 


function  [ f e, ee, f eu, zhc]  =  Felastic ( z,  zeq,  R, est , H, seg, gama ) 
f e ( 1 : length ( z) ) =0; 
ee ( 1 : length ( z) )  =0; 

K  =  . 545+ . 41* .25;  %  related  to  poisson  ration  v=.3 
KK=4/3*est/ ( 1- . 25A2 ) /pi A2; 

EE=est/ (1-.25A2)  ; 
zhc= (pi*K*H/2/EE) A  2  *  R  ; 
for  k=l : length ( z) 
zh=z ( k) -zeq; 
if  zh  >  0 

if  zh  <=  zhc 

fe  (k)  =4/3*EE*sqrt  (R)  *zh/'l  .5; 
feu (k) =2*R*zh*pi*H; 
ee(k)=.5*zh* (fe(k)-feu(k) ) ; 

elseif  zh  >  zhc  &  zh  <=  6*zhc 

fe (k) =1 . 03* (zh/zhc) A1 . 425* 4/3*EE*sqrt (R) *zhcA1.5; 

%  ee(k)=l/2*fe(k)*zh; 

feu (k) =2*R*zh*pi*H; 
ee(k)=.5*zh* (fe(k)-feu(k) ) ; 

elseif  zh  >  6*zhc  &  zh  <=  110*zhc 

%  fe(k)=2*pi*zh*R*H; 

f e ( k) =1 . 4* (zh/zhc) A1 ,263*4/3*EE*sqrt (R) *zhcAl .5; 

% 

ee (k) =2/5* (4/3*EE*sqrt (R) *zhcAl .5) A (5/3) /KKA (2/3) /RA (1/3) ; 

I  ee(k)=l/2*fe(k)*zh; 

feu(k)=2*R*zh*pi*H; 
ee(k)=.5*zh* (fe(k)-feu(k) ) ; 

elseif  zh  >  110*zhc  &  zh  <=  R 

fe (k)=3/K* (zh/zhc) Al*4/3*EE*sqrt (R) *zhcA1.5; 
f eu ( k) =2*R*zh*pi*H; 
ee ( k)  = . 5*  zh* (fe(k) -feu(k) ) ; 

%  ee(k)=l/2*fe(k) *zh-l/2*fe (k) * (zh-110*zhc) ; 

% 

ee (k) =2/5* (4/3*EE*sqrt (R) *zhcA1.5) A (5/3) /KKA (2/3) /RA (1/3) ; 
elseif  zh  >  110*zhc  &  zh  >=  R 
f eu ( k) =2*R* zh*pi*H; 

fe(k)=3/K* (zh/zhc) Al*4/3*EE*sqrt (R) *zhcAl .5; 
ee(k)=.5*zh*(fe(k)-feu(k)); 

%  fe (k) =pi*RA2*H; 

%  ee ( k) =pi*H* (RA3/3-l/3* (R-zh) A3-RA2*zh) ; 

%  ee (k) =1/2* fe (k) *zh-l/2*fe (k) * ( zh-110* zhc) ; 
end 

end 

end 


function  ( f , e] =fcapx ( z, ps2 , x, zeq, R, gama ,thl,th2,RH,Vl,Rg,T,seg, zhc, dz ) 
W(1 : length (x) ) =0; 


fh  ( 1 : length (x) ) =0 ; 
e  ( 1 : length (x) ) =0 ; 
f  ( 1 : length (x) ) =0; 

%  Finding  the  capillary  force 


and  surface  energy  as  a  function  of  x 


o.  a  Q.  o.  ' 
o  %  ?  ~6  ' 

~o  -o  ^  x>  o'o'o'o'e  O  f  fl  f  t  o  o  o  o  o  •&  "S  •p  o'&'&'S 

for  k=l : length (x) 


h=zeq-x ( k) ; 
if  h  <=  -R 
h=-R; 

end 

phi=ps2 ( k) ; 

%%%% 

%%%%%%%%%%%%%%%%%%% 

rl=( ( (R* (1-cos (phi) )+h) / (cos ( thl+phi ) +cos ( th2) ) ) ) ; 
if  phi  ~=  0 

fh(k)=2*pi*R*sin(phi)*sin(phi+thl)-(- 
1* (cos (th2) +cos (thl+phi) ) / (R- 

R*cos (phi ) +h) +sin (thl+phi )/sin(phi)/R)*pi*(R*sin (phi) ) ~2 ; 

Esls=-1  * pi* cos (th2) * (R*sin(phi)+rl*sin (thl+phi) - 
rl*sin (th2) ) A2; 

Esla=-l*2*pi*R''2*  (1-cos  (phi)  )  *cos  (thl )  ; 
if  h  <  0 

c=sqrt (R"2- (R+h) ^2) ; 

Esls=Esls+pi*c/'2*cos  (th2)  ; 

Esla=Esla+pi*2*R*h*-l*cos (thl ) ; 

fh (k) =2*pi*R*sin (phi ) *sin (phi+thl ) - (- 
1* (cos (th2) +cos (thl+phi) ) / (R- 

R*cos  (phi  )  +h)  +sin  (thl+phi )  /  sin  (phi )  /R)  *pi*  (  (R*sin  (phi  )  )/'2-c'~2)  ; 
end 

Elv=2*pi* ( (cos ( th2 ) +cos (thl+phi ) ) * (sin ( phi ) +sin (thl+phi) ) + (thl+phi+th2) 
/2-l/4*sin (2*thl+2*phi )-l/4*sin(2*th2) ) *rl^2; 

W (k) =Elv+Esls+Esla; 
if  h  <=  -R 

W (k) =W (k-1) ; 
fh (k) =fh ( k— 1 ) ; 

end 

end 

end 


f=fh;  %  force  density  as  function  of  z 
e=W;  %  energy  density  as  function  of  z 


function  [ fv, ev) =fvdw (Ahl , Ahg, z,zeq,R,x,ps2,thl, th2 , seg, gama, zhc, dz) 
do= . 17e-9/seg; 
fv ( 1 : length ( z) ) =0; 
ev ( 1 : length ( z) ) =0; 
for  i=l : length ( z ) 
h=zeq-z (i ) ; 
if  h  <=  -R 


i  =  <  3 


h=-R; 

end 

D=h; 

phi=ps2  ( i  ) .; 

c.Q.Q,o 
"6  'o  "o  "o 

rl= ( ( ( (R*  (1-cos (phi) ) +h) / (cos (thl+phi) +cos (th2) )  )  )  )  ; 
zcc=rl*  (cos (phi+thl) +cos (th2) ) -D; 
if  rl  ~=  0  &  phi  ~=  0  &  D  ~=  0 

if  D  >=  do 

o 

"c 

fv (i) =R* (Ahg- 

Ahl)  / 6/ ( zcc+D) A2* ( (3*zcc+D) / (zcc+D) +Ahl* (zcc+D) A2/ ( Ahg-Ahl ) /DA2 ) ; 
ev (i) =-l* (Ahg- 

Ahl  )  *  R/6/ ( zcc+D) * ( (2*zcc+D)/( zcc+D) +Ahl* ( zcc+D) / (Ahl-Ahg) / D )  ; 
'elseif  D  <  do 

fv (i) =R* (Ahg- 

Ahl)  /  6/  (zcc+D)  A2* ( (3* zcc+D) / (zcc+D) +Ahl * (zcc+D) A2* (3*do-2*D) / (Ahg- 
Ahl)  /DA2)  ; 

ev ( i ) =-l * (Ahg-Ahl )*R/6/(zcc+D)* ( (2*zcc+D)/ (zcc+D)- 
Ahl* (zcc+D) * (2*do+D) / (Ahl-Ahg) /doA2) ; 
else 

ev  ( i ) =0; 
fv (i) =0; 


end 

end 

if  phi  ==  0  &  h  >=  0 

fv ( i ) =Ahg/6* (l/hA2) *R; 
ev (i) =-Ahg/6* (1/h) *R; 

end 

if  phi  ==  0  &  h  <  do 

f v ( i ) — Ahg / 6  * (1/ (do) A2) *R; 
ev (i) =-Ahg/6* (1/ (do) ) *R; 

end 

end 

fv=abs ( fv) ; 

%  plot ( fv) 

%  pause 


function  A=integ(F,x) 
sum=0; 

for  i=2 : length (x) 

sum=sunt+  (F(i)  +F(i  —  1)  )  / 2*  (x(i)  -x(i-l)  )  ; 

end 

A=sum; 


